首页 > 代码库 > temp、
temp、
1 #include <map> 2 #include <set> 3 #include <cmath> 4 #include <queue> 5 #include <stack> 6 #include <cstdio> 7 #include <string> 8 #include <vector> 9 #include <cstdlib> 10 #include <cstring> 11 #include <sstream> 12 #include <iostream> 13 #include <algorithm> 14 #include <functional> 15 using namespace std; 16 #define rep(i,a,n) for (int i=a;i<n;i++) 17 #define per(i,a,n) for (int i=n-1;i>=a;i--) 18 #define all(x) (x).begin(),(x).end() 19 #define pb push_back 20 #define mp make_pair 21 #define lson l,m,rt<<1 22 #define rson m+1,r,rt<<1|1 23 typedef long long ll; 24 typedef vector<int> VI; 25 typedef pair<int, int> PII; 26 const ll MOD = 1e9 + 7; 27 const int INF = 0x3f3f3f3f; 28 const int MAXN = 1e5 + 7; 29 // head 30 /* 31 计算几何 32 33 目录 34 ㈠ 点的基本运算 35 1. 平面上两点之间距离 1 36 2. 判断两点是否重合 1 37 3. 矢量叉乘 1 38 4. 矢量点乘 2 39 5. 判断点是否在线段上 2 40 6. 求一点饶某点旋转后的坐标 2 41 7. 求矢量夹角 2 42 43 ㈡ 线段及直线的基本运算 44 1. 点与线段的关系 3 45 2. 求点到线段所在直线垂线的垂足 4 46 3. 点到线段的最近点 4 47 4. 点到线段所在直线的距离 4 48 5. 点到折线集的最近距离 4 49 6. 判断圆是否在多边形内 5 50 7. 求矢量夹角余弦 5 51 8. 求线段之间的夹角 5 52 9. 判断线段是否相交 6 53 10.判断线段是否相交但不交在端点处 6 54 11.求线段所在直线的方程 6 55 12.求直线的斜率 7 56 13.求直线的倾斜角 7 57 14.求点关于某直线的对称点 7 58 15.判断两条直线是否相交及求直线交点 7 59 16.判断线段是否相交，如果相交返回交点 7 60 61 ㈢ 多边形常用算法模块 62 1. 判断多边形是否简单多边形 8 63 2. 检查多边形顶点的凸凹性 9 64 3. 判断多边形是否凸多边形 9 65 4. 求多边形面积 9 66 5. 判断多边形顶点的排列方向，方法一 10 67 6. 判断多边形顶点的排列方向，方法二 10 68 7. 射线法判断点是否在多边形内 10 69 8. 判断点是否在凸多边形内 11 70 9. 寻找点集的graham算法 12 71 10.寻找点集凸包的卷包裹法 13 72 11.判断线段是否在多边形内 14 73 12.求简单多边形的重心 15 74 13.求凸多边形的重心 17 75 14.求肯定在给定多边形内的一个点 17 76 15.求从多边形外一点出发到该多边形的切线 18 77 16.判断多边形的核是否存在 19 78 79 ㈣ 圆的基本运算 80 .点是否在圆内 20 81 .求不共线的三点所确定的圆 21 82 83 ㈤ 矩形的基本运算 84 1.已知矩形三点坐标，求第4点坐标 22 85 86 ㈥ 常用算法的描述 22 87 88 ㈦ 补充 89 1．两圆关系： 24 90 2．判断圆是否在矩形内： 24 91 3．点到平面的距离： 25 92 4．点是否在直线同侧： 25 93 5．镜面反射线： 25 94 6．矩形包含： 26 95 7．两圆交点： 27 96 8．两圆公共面积： 28 97 9. 圆和直线关系： 29 98 10. 内切圆： 30 99 11. 求切点： 31 100 12. 线段的左右旋： 31 101 13．公式： 32 102 * / 103 /* 需要包含的头文件 */ 104 105 #include <cmath > 106 107 /* 常用的常量定义 */ 108 const double INF = 1e200; 109 const double EP = 1e-10; 110 const int MAXV = 300; 111 const double PI = 3.14159265; 112 113 /* 基本几何结构 */ 114 struct POINT 115 { 116 double x; 117 double y; 118 POINT(double a = 0, double b = 0) { x = a; y = b; } //constructor 119 }; 120 struct LINESEG 121 { 122 POINT s; 123 POINT e; 124 LINESEG(POINT a, POINT b) { s = a; e = b; } 125 LINESEG() { } 126 }; 127 struct LINE // 直线的解析方程 a*x+b*y+c=0 为统一表示，约定 a >= 0 128 { 129 double a; 130 double b; 131 double c; 132 LINE(double d1 = 1, double d2 = -1, double d3 = 0) { a = d1; b = d2; c = d3; } 133 }; 134 135 /********************** 136 * * 137 * 点的基本运算 * 138 * * 139 **********************/ 140 141 double dist(POINT p1, POINT p2) // 返回两点之间欧氏距离 142 { 143 return(sqrt((p1.x - p2.x)*(p1.x - p2.x) + (p1.y - p2.y)*(p1.y - p2.y))); 144 } 145 bool equal_point(POINT p1, POINT p2) // 判断两个点是否重合 146 { 147 return ((abs(p1.x - p2.x)<EP) && (abs(p1.y - p2.y)<EP)); 148 } 149 /****************************************************************************** 150 r=multiply(sp,ep,op),得到(sp-op)和(ep-op)的叉积 151 r>0：ep在矢量opsp的逆时针方向； 152 r=0：opspep三点共线； 153 r<0：ep在矢量opsp的顺时针方向 154 *******************************************************************************/ 155 double multiply(POINT sp, POINT ep, POINT op) 156 { 157 return((sp.x - op.x)*(ep.y - op.y) - (ep.x - op.x)*(sp.y - op.y)); 158 } 159 /* 160 r=dotmultiply(p1,p2,op),得到矢量(p1-op)和(p2-op)的点积，如果两个矢量都非零矢量 161 r<0：两矢量夹角为钝角； 162 r=0：两矢量夹角为直角； 163 r>0：两矢量夹角为锐角 164 *******************************************************************************/ 165 double dotmultiply(POINT p1, POINT p2, POINT p0) 166 { 167 return ((p1.x - p0.x)*(p2.x - p0.x) + (p1.y - p0.y)*(p2.y - p0.y)); 168 } 169 /****************************************************************************** 170 判断点p是否在线段l上 171 条件：(p在线段l所在的直线上) && (点p在以线段l为对角线的矩形内) 172 *******************************************************************************/ 173 bool online(LINESEG l, POINT p) 174 { 175 return((multiply(l.e, p, l.s) == 0) && (((p.x - l.s.x)*(p.x - l.e.x) <= 0) && ((p.y - l.s.y)*(p.y - l.e.y) <= 0))); 176 } 177 // 返回点p以点o为圆心逆时针旋转alpha(单位：弧度)后所在的位置 178 POINT rotate(POINT o, double alpha, POINT p) 179 { 180 POINT tp; 181 p.x -= o.x; 182 p.y -= o.y; 183 tp.x = p.x*cos(alpha) - p.y*sin(alpha) + o.x; 184 tp.y = p.y*cos(alpha) + p.x*sin(alpha) + o.y; 185 return tp; 186 } 187 /* 返回顶角在o点，起始边为os，终止边为oe的夹角(单位：弧度) 188 角度小于pi，返回正值 189 角度大于pi，返回负值 190 可以用于求线段之间的夹角 191 原理： 192 r = dotmultiply(s,e,o) / (dist(o,s)*dist(o,e)) 193 r‘= multiply(s,e,o) 194 195 r >= 1 angle = 0; 196 r <= -1 angle = -PI 197 -1<r<1 && r‘>0 angle = arccos(r) 198 -1<r<1 && r‘<=0 angle = -arccos(r) 199 */ 200 double angle(POINT o, POINT s, POINT e) 201 { 202 double cosfi, fi, norm; 203 double dsx = s.x - o.x; 204 double dsy = s.y - o.y; 205 double dex = e.x - o.x; 206 double dey = e.y - o.y; 207 208 cosfi = dsx*dex + dsy*dey; 209 norm = (dsx*dsx + dsy*dsy)*(dex*dex + dey*dey); 210 cosfi /= sqrt(norm); 211 212 if (cosfi >= 1.0) return 0; 213 if (cosfi <= -1.0) return -3.1415926; 214 215 fi = acos(cosfi); 216 if (dsx*dey - dsy*dex>0) return fi; // 说明矢量os 在矢量 oe的顺时针方向 217 return -fi; 218 } 219 /***************************** 220 * * 221 * 线段及直线的基本运算 * 222 * * 223 \*****************************/ 224 225 /* 判断点与线段的关系,用途很广泛 226 本函数是根据下面的公式写的，P是点C到线段AB所在直线的垂足 227 228 AC dot AB 229 r = --------- 230 ||AB||^2 231 (Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay) 232 = ------------------------------- 233 L^2 234 235 r has the following meaning: 236 237 r=0 P = A 238 r=1 P = B 239 r<0 P is on the backward extension of AB 240 r>1 P is on the forward extension of AB 241 0<r<1 P is interior to AB 242 */ 243 double relation(POINT p, LINESEG l) 244 { 245 LINESEG tl; 246 tl.s = l.s; 247 tl.e = p; 248 return dotmultiply(tl.e, l.e, l.s) / (dist(l.s, l.e)*dist(l.s, l.e)); 249 } 250 // 求点C到线段AB所在直线的垂足 P 251 POINT perpendicular(POINT p, LINESEG l) 252 { 253 double r = relation(p, l); 254 POINT tp; 255 tp.x = l.s.x + r*(l.e.x - l.s.x); 256 tp.y = l.s.y + r*(l.e.y - l.s.y); 257 return tp; 258 } 259 /* 求点p到线段l的最短距离,并返回线段上距该点最近的点np 260 注意：np是线段l上到点p最近的点，不一定是垂足 */ 261 double ptolinesegdist(POINT p, LINESEG l, POINT &np) 262 { 263 double r = relation(p, l); 264 if (r<0) 265 { 266 np = l.s; 267 return dist(p, l.s); 268 } 269 if (r>1) 270 { 271 np = l.e; 272 return dist(p, l.e); 273 } 274 np = perpendicular(p, l); 275 return dist(p, np); 276 } 277 // 求点p到线段l所在直线的距离,请注意本函数与上个函数的区别 278 double ptoldist(POINT p, LINESEG l) 279 { 280 return abs(multiply(p, l.e, l.s)) / dist(l.s, l.e); 281 } 282 /* 计算点到折线集的最近距离,并返回最近点. 283 注意：调用的是ptolineseg()函数 */ 284 double ptopointset(int vcount, POINT pointset[], POINT p, POINT &q) 285 { 286 int i; 287 double cd = double(INF), td; 288 LINESEG l; 289 POINT tq, cq; 290 291 for (i = 0; i<vcount - 1; i++) 292 { 293 l.s = pointset[i]; 294 295 l.e = pointset[i + 1]; 296 td = ptolinesegdist(p, l, tq); 297 if (td<cd) 298 { 299 cd = td; 300 cq = tq; 301 } 302 } 303 q = cq; 304 return cd; 305 } 306 /* 判断圆是否在多边形内.ptolineseg()函数的应用2 */ 307 bool CircleInsidePolygon(int vcount, POINT center, double radius, POINT polygon[]) 308 { 309 POINT q; 310 double d; 311 q.x = 0; 312 q.y = 0; 313 d = ptopointset(vcount, polygon, center, q); 314 if (d<radius || fabs(d - radius)<EP) 315 return true; 316 else 317 return false; 318 } 319 /* 返回两个矢量l1和l2的夹角的余弦(-1 --- 1)注意：如果想从余弦求夹角的话，注意反余弦函数的定义域是从 0到pi */ 320 double cosine(LINESEG l1, LINESEG l2) 321 { 322 return (((l1.e.x - l1.s.x)*(l2.e.x - l2.s.x) + 323 (l1.e.y - l1.s.y)*(l2.e.y - l2.s.y)) / (dist(l1.e, l1.s)*dist(l2.e, l2.s))); 324 } 325 // 返回线段l1与l2之间的夹角 单位：弧度 范围(-pi，pi) 326 double lsangle(LINESEG l1, LINESEG l2) 327 { 328 POINT o, s, e; 329 o.x = o.y = 0; 330 s.x = l1.e.x - l1.s.x; 331 s.y = l1.e.y - l1.s.y; 332 e.x = l2.e.x - l2.s.x; 333 e.y = l2.e.y - l2.s.y; 334 return angle(o, s, e); 335 } 336 // 如果线段u和v相交(包括相交在端点处)时，返回true 337 // 338 //判断P1P2跨立Q1Q2的依据是：( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) >= 0。 339 //判断Q1Q2跨立P1P2的依据是：( Q1 - P1 ) × ( P2 - P1 ) * ( P2 - P1 ) × ( Q2 - P1 ) >= 0。 340 bool intersect(LINESEG u, LINESEG v) 341 { 342 return((max(u.s.x, u.e.x) >= min(v.s.x, v.e.x)) && //排斥实验 343 (max(v.s.x, v.e.x) >= min(u.s.x, u.e.x)) && 344 (max(u.s.y, u.e.y) >= min(v.s.y, v.e.y)) && 345 (max(v.s.y, v.e.y) >= min(u.s.y, u.e.y)) && 346 (multiply(v.s, u.e, u.s)*multiply(u.e, v.e, u.s) >= 0) && //跨立实验 347 (multiply(u.s, v.e, v.s)*multiply(v.e, u.e, v.s) >= 0)); 348 } 349 // (线段u和v相交)&&(交点不是双方的端点) 时返回true 350 bool intersect_A(LINESEG u, LINESEG v) 351 { 352 return ((intersect(u, v)) && 353 (!online(u, v.s)) && 354 (!online(u, v.e)) && 355 (!online(v, u.e)) && 356 (!online(v, u.s))); 357 } 358 // 线段v所在直线与线段u相交时返回true；方法：判断线段u是否跨立线段v 359 bool intersect_l(LINESEG u, LINESEG v) 360 { 361 return multiply(u.s, v.e, v.s)*multiply(v.e, u.e, v.s) >= 0; 362 } 363 // 根据已知两点坐标，求过这两点的直线解析方程： a*x+b*y+c = 0 (a >= 0) 364 LINE makeline(POINT p1, POINT p2) 365 { 366 LINE tl; 367 int sign = 1; 368 tl.a = p2.y - p1.y; 369 if (tl.a<0) 370 { 371 sign = -1; 372 tl.a = sign*tl.a; 373 } 374 tl.b = sign*(p1.x - p2.x); 375 tl.c = sign*(p1.y*p2.x - p1.x*p2.y); 376 return tl; 377 } 378 // 根据直线解析方程返回直线的斜率k,水平线返回 0,竖直线返回 1e200 379 double slope(LINE l) 380 { 381 if (abs(l.a) < 1e-20) 382 return 0; 383 if (abs(l.b) < 1e-20) 384 return INF; 385 return -(l.a / l.b); 386 } 387 // 返回直线的倾斜角alpha ( 0 - pi) 388 double alpha(LINE l) 389 { 390 if (abs(l.a)< EP) 391 return 0; 392 if (abs(l.b)< EP) 393 return PI / 2; 394 double k = slope(l); 395 if (k>0) 396 return atan(k); 397 else 398 return PI + atan(k); 399 } 400 // 求点p关于直线l的对称点 401 POINT symmetry(LINE l, POINT p) 402 { 403 POINT tp; 404 tp.x = ((l.b*l.b - l.a*l.a)*p.x - 2 * l.a*l.b*p.y - 2 * l.a*l.c) / (l.a*l.a + l.b*l.b); 405 tp.y = ((l.a*l.a - l.b*l.b)*p.y - 2 * l.a*l.b*p.x - 2 * l.b*l.c) / (l.a*l.a + l.b*l.b); 406 return tp; 407 } 408 // 如果两条直线 l1(a1*x+b1*y+c1 = 0), l2(a2*x+b2*y+c2 = 0)相交，返回true，且返回交点p 409 bool lineintersect(LINE l1, LINE l2, POINT &p) // 是 L1，L2 410 { 411 double d = l1.a*l2.b - l2.a*l1.b; 412 if (abs(d)<EP) // 不相交 413 return false; 414 p.x = (l2.c*l1.b - l1.c*l2.b) / d; 415 p.y = (l2.a*l1.c - l1.a*l2.c) / d; 416 return true; 417 } 418 // 如果线段l1和l2相交，返回true且交点由(inter)返回，否则返回false 419 bool intersection(LINESEG l1, LINESEG l2, POINT &inter) 420 { 421 LINE ll1, ll2; 422 ll1 = makeline(l1.s, l1.e); 423 ll2 = makeline(l2.s, l2.e); 424 if (lineintersect(ll1, ll2, inter)) 425 return online(l1, inter) && online(l2, inter); 426 else 427 return false; 428 } 429 430 /****************************** 431 * * 432 * 多边形常用算法模块 * 433 * * 434 \******************************/ 435 436 // 如果无特别说明，输入多边形顶点要求按逆时针排列 437 438 /* 439 返回值：输入的多边形是简单多边形，返回true 440 要 求：输入顶点序列按逆时针排序 441 说 明：简单多边形定义： 442 1：循环排序中相邻线段对的交是他们之间共有的单个点 443 2：不相邻的线段不相交 444 本程序默认第一个条件已经满足 445 */ 446 bool issimple(int vcount, POINT polygon[]) 447 { 448 int i, cn; 449 LINESEG l1, l2; 450 for (i = 0; i<vcount; i++) 451 { 452 l1.s = polygon[i]; 453 l1.e = polygon[(i + 1) % vcount]; 454 cn = vcount - 3; 455 while (cn) 456 { 457 l2.s = polygon[(i + 2) % vcount]; 458 l2.e = polygon[(i + 3) % vcount]; 459 if (intersect(l1, l2)) 460 break; 461 cn--; 462 } 463 if (cn) 464 return false; 465 } 466 return true; 467 } 468 // 返回值：按输入顺序返回多边形顶点的凸凹性判断，bc[i]=1,iff:第i个顶点是凸顶点 469 void checkconvex(int vcount, POINT polygon[], bool bc[]) 470 { 471 int i, index = 0; 472 POINT tp = polygon[0]; 473 for (i = 1; i<vcount; i++) // 寻找第一个凸顶点 474 { 475 if (polygon[i].y<tp.y || (polygon[i].y == tp.y&&polygon[i].x<tp.x)) 476 { 477 tp = polygon[i]; 478 index = i; 479 } 480 } 481 int count = vcount - 1; 482 bc[index] = 1; 483 while (count) // 判断凸凹性 484 { 485 if (multiply(polygon[(index + 1) % vcount], polygon[(index + 2) % vcount], polygon[index]) >= 0) 486 bc[(index + 1) % vcount] = 1; 487 else 488 bc[(index + 1) % vcount] = 0; 489 index++; 490 count--; 491 } 492 } 493 // 返回值：多边形polygon是凸多边形时，返回true 494 bool isconvex(int vcount, POINT polygon[]) 495 { 496 bool bc[MAXV]; 497 checkconvex(vcount, polygon, bc); 498 for (int i = 0; i<vcount; i++) // 逐一检查顶点，是否全部是凸顶点 499 if (!bc[i]) 500 return false; 501 return true; 502 } 503 // 返回多边形面积(signed)；输入顶点按逆时针排列时，返回正值；否则返回负值 504 double area_of_polygon(int vcount, POINT polygon[]) 505 { 506 int i; 507 double s; 508 if (vcount<3) 509 return 0; 510 s = polygon[0].y*(polygon[vcount - 1].x - polygon[1].x); 511 for (i = 1; i<vcount; i++) 512 s += polygon[i].y*(polygon[(i - 1)].x - polygon[(i + 1) % vcount].x); 513 return s / 2; 514 } 515 // 如果输入顶点按逆时针排列，返回true 516 bool isconterclock(int vcount, POINT polygon[]) 517 { 518 return area_of_polygon(vcount, polygon)>0; 519 } 520 // 另一种判断多边形顶点排列方向的方法 521 bool isccwize(int vcount, POINT polygon[]) 522 { 523 int i, index; 524 POINT a, b, v; 525 v = polygon[0]; 526 index = 0; 527 for (i = 1; i<vcount; i++) // 找到最低且最左顶点，肯定是凸顶点 528 { 529 if (polygon[i].y<v.y || polygon[i].y == v.y && polygon[i].x<v.x) 530 { 531 index = i; 532 } 533 } 534 a = polygon[(index - 1 + vcount) % vcount]; // 顶点v的前一顶点 535 b = polygon[(index + 1) % vcount]; // 顶点v的后一顶点 536 return multiply(v, b, a)>0; 537 } 538 /******************************************************************************************** 539 射线法判断点q与多边形polygon的位置关系，要求polygon为简单多边形，顶点逆时针排列 540 如果点在多边形内： 返回0 541 如果点在多边形边上： 返回1 542 如果点在多边形外： 返回2 543 *********************************************************************************************/ 544 int insidepolygon(int vcount, POINT Polygon[], POINT q) 545 { 546 int c = 0, i, n; 547 LINESEG l1, l2; 548 bool bintersect_a, bonline1, bonline2, bonline3; 549 double r1, r2; 550 551 l1.s = q; 552 l1.e = q; 553 l1.e.x = double(INF); 554 n = vcount; 555 for (i = 0; i<vcount; i++) 556 { 557 l2.s = Polygon[i]; 558 l2.e = Polygon[(i + 1) % n]; 559 if (online(l2, q)) 560 return 1; // 如果点在边上，返回1 561 if ((bintersect_a = intersect_A(l1, l2)) || // 相交且不在端点 562 ((bonline1 = online(l1, Polygon[(i + 1) % n])) && // 第二个端点在射线上 563 ((!(bonline2 = online(l1, Polygon[(i + 2) % n]))) && /* 前一个端点和后一个端点在射线两侧 */ 564 ((r1 = multiply(Polygon[i], Polygon[(i + 1) % n], l1.s)*multiply(Polygon[(i + 1) % n], Polygon[(i + 2) % n], l1.s))>0) || 565 (bonline3 = online(l1, Polygon[(i + 2) % n])) && /* 下一条边是水平线，前一个端点和后一个端点在射线两侧 */ 566 ((r2 = multiply(Polygon[i], Polygon[(i + 2) % n], l1.s)*multiply(Polygon[(i + 2) % n], 567 Polygon[(i + 3) % n], l1.s))>0) 568 ) 569 ) 570 ) c++; 571 } 572 if (c % 2 == 1) 573 return 0; 574 else 575 return 2; 576 } 577 //点q是凸多边形polygon内时，返回true；注意：多边形polygon一定要是凸多边形 578 bool InsideConvexPolygon(int vcount, POINT polygon[], POINT q) // 可用于三角形！ 579 { 580 POINT p; 581 LINESEG l; 582 int i; 583 p.x = 0; p.y = 0; 584 for (i = 0; i<vcount; i++) // 寻找一个肯定在多边形polygon内的点p：多边形顶点平均值 585 { 586 p.x += polygon[i].x; 587 p.y += polygon[i].y; 588 } 589 p.x /= vcount; 590 p.y /= vcount; 591 592 for (i = 0; i<vcount; i++) 593 { 594 l.s = polygon[i]; l.e = polygon[(i + 1) % vcount]; 595 if (multiply(p, l.e, l.s)*multiply(q, l.e, l.s)<0) /* 点p和点q在边l的两侧，说明点q肯定在多边形外 */ 596 break; 597 } 598 return (i == vcount); 599 } 600 /********************************************** 601 寻找凸包的graham 扫描法 602 PointSet为输入的点集； 603 ch为输出的凸包上的点集，按照逆时针方向排列; 604 n为PointSet中的点的数目 605 len为输出的凸包上的点的个数 606 **********************************************/ 607 void Graham_scan(POINT PointSet[], POINT ch[], int n, int &len) 608 { 609 int i, j, k = 0, top = 2; 610 POINT tmp; 611 // 选取PointSet中y坐标最小的点PointSet[k]，如果这样的点有多个，则取最左边的一个 612 for (i = 1; i<n; i++) 613 if (PointSet[i].y<PointSet[k].y || (PointSet[i].y == PointSet[k].y) && (PointSet[i].x<PointSet[k].x)) 614 k = i; 615 tmp = PointSet[0]; 616 PointSet[0] = PointSet[k]; 617 PointSet[k] = tmp; // 现在PointSet中y坐标最小的点在PointSet[0] 618 for (i = 1; i<n - 1; i++) /* 对顶点按照相对PointSet[0]的极角从小到大进行排序，极角相同的按照距离PointSet[0]从近到远进行排序 */ 619 { 620 k = i; 621 for (j = i + 1; j<n; j++) 622 if (multiply(PointSet[j], PointSet[k], PointSet[0])>0 || // 极角更小 623 (multiply(PointSet[j], PointSet[k], PointSet[0]) == 0) && /* 极角相等，距离更短 */ 624 dist(PointSet[0], PointSet[j])<dist(PointSet[0], PointSet[k]) 625 ) 626 k = j; 627 tmp = PointSet[i]; 628 PointSet[i] = PointSet[k]; 629 PointSet[k] = tmp; 630 } 631 ch[0] = PointSet[0]; 632 ch[1] = PointSet[1]; 633 ch[2] = PointSet[2]; 634 for (i = 3; i<n; i++) 635 { 636 while (multiply(PointSet[i], ch[top], ch[top - 1]) >= 0) 637 top--; 638 ch[++top] = PointSet[i]; 639 } 640 len = top + 1; 641 } 642 // 卷包裹法求点集凸壳，参数说明同graham算法 643 void ConvexClosure(POINT PointSet[], POINT ch[], int n, int &len) 644 { 645 int top = 0, i, index, first; 646 double curmax, curcos, curdis; 647 POINT tmp; 648 LINESEG l1, l2; 649 bool use[MAXV]; 650 tmp = PointSet[0]; 651 index = 0; 652 // 选取y最小点，如果多于一个，则选取最左点 653 for (i = 1; i<n; i++) 654 { 655 if (PointSet[i].y<tmp.y || PointSet[i].y == tmp.y&&PointSet[i].x<tmp.x) 656 { 657 index = i; 658 } 659 use[i] = false; 660 } 661 tmp = PointSet[index]; 662 first = index; 663 use[index] = true; 664 665 index = -1; 666 ch[top++] = tmp; 667 tmp.x -= 100; 668 l1.s = tmp; 669 l1.e = ch[0]; 670 l2.s = ch[0]; 671 672 while (index != first) 673 { 674 curmax = -100; 675 curdis = 0; 676 // 选取与最后一条确定边夹角最小的点，即余弦值最大者 677 for (i = 0; i<n; i++) 678 { 679 if (use[i])continue; 680 l2.e = PointSet[i]; 681 curcos = cosine(l1, l2); // 根据cos值求夹角余弦，范围在 （-1 -- 1 ） 682 if (curcos>curmax || fabs(curcos - curmax)<1e-6 && dist(l2.s, l2.e)>curdis) 683 { 684 curmax = curcos; 685 index = i; 686 curdis = dist(l2.s, l2.e); 687 } 688 } 689 use[first] = false; //清空第first个顶点标志，使最后能形成封闭的hull 690 use[index] = true; 691 ch[top++] = PointSet[index]; 692 l1.s = ch[top - 2]; 693 l1.e = ch[top - 1]; 694 l2.s = ch[top - 1]; 695 } 696 len = top - 1; 697 } 698 /********************************************************************************************* 699 判断线段是否在简单多边形内(注意：如果多边形是凸多边形，下面的算法可以化简) 700 必要条件一：线段的两个端点都在多边形内； 701 必要条件二：线段和多边形的所有边都不内交； 702 用途： 1. 判断折线是否在简单多边形内 703 2. 判断简单多边形是否在另一个简单多边形内 704 **********************************************************************************************/ 705 bool LinesegInsidePolygon(int vcount, POINT polygon[], LINESEG l) 706 { 707 // 判断线端l的端点是否不都在多边形内 708 if (!insidepolygon(vcount, polygon, l.s) || !insidepolygon(vcount, polygon, l.e)) 709 return false; 710 int top = 0, i, j; 711 POINT PointSet[MAXV], tmp; 712 LINESEG s; 713 714 for (i = 0; i<vcount; i++) 715 { 716 s.s = polygon[i]; 717 s.e = polygon[(i + 1) % vcount]; 718 if (online(s, l.s)) //线段l的起始端点在线段s上 719 PointSet[top++] = l.s; 720 else if (online(s, l.e)) //线段l的终止端点在线段s上 721 PointSet[top++] = l.e; 722 else 723 { 724 if (online(l, s.s)) //线段s的起始端点在线段l上 725 PointSet[top++] = s.s; 726 else if (online(l, s.e)) // 线段s的终止端点在线段l上 727 PointSet[top++] = s.e; 728 else 729 { 730 if (intersect(l, s)) // 这个时候如果相交，肯定是内交，返回false 731 return false; 732 } 733 } 734 } 735 736 for (i = 0; i<top - 1; i++) /* 冒泡排序，x坐标小的排在前面；x坐标相同者，y坐标小的排在前面 */ 737 { 738 for (j = i + 1; j<top; j++) 739 { 740 if (PointSet[i].x>PointSet[j].x || fabs(PointSet[i].x - PointSet[j].x)<EP && PointSet[i].y>PointSet[j].y) 741 { 742 tmp = PointSet[i]; 743 PointSet[i] = PointSet[j]; 744 PointSet[j] = tmp; 745 } 746 } 747 } 748 749 for (i = 0; i<top - 1; i++) 750 { 751 tmp.x = (PointSet[i].x + PointSet[i + 1].x) / 2; //得到两个相邻交点的中点 752 tmp.y = (PointSet[i].y + PointSet[i + 1].y) / 2; 753 if (!insidepolygon(vcount, polygon, tmp)) 754 return false; 755 } 756 return true; 757 } 758 /********************************************************************************************* 759 求任意简单多边形polygon的重心 760 需要调用下面几个函数： 761 void AddPosPart(); 增加右边区域的面积 762 void AddNegPart(); 增加左边区域的面积 763 void AddRegion(); 增加区域面积 764 在使用该程序时，如果把xtr,ytr,wtr,xtl,ytl,wtl设成全局变量就可以使这些函数的形式得到化简, 765 但要注意函数的声明和调用要做相应变化 766 **********************************************************************************************/ 767 void AddPosPart(double x, double y, double w, double &xtr, double &ytr, double &wtr) 768 { 769 if (abs(wtr + w)<1e-10) return; // detect zero regions 770 xtr = (wtr*xtr + w*x) / (wtr + w); 771 ytr = (wtr*ytr + w*y) / (wtr + w); 772 wtr = w + wtr; 773 return; 774 } 775 void AddNegPart(double x, double y, double w, double &xtl, double &ytl, double &wtl) 776 { 777 if (abs(wtl + w)<1e-10) 778 return; // detect zero regions 779 780 xtl = (wtl*xtl + w*x) / (wtl + w); 781 ytl = (wtl*ytl + w*y) / (wtl + w); 782 wtl = w + wtl; 783 return; 784 } 785 void AddRegion(double x1, double y1, double x2, double y2, double &xtr, double &ytr, 786 double &wtr, double &xtl, double &ytl, double &wtl) 787 { 788 if (abs(x1 - x2)< 1e-10) 789 return; 790 791 if (x2 > x1) 792 { 793 AddPosPart((x2 + x1) / 2, y1 / 2, (x2 - x1) * y1, xtr, ytr, wtr); /* rectangle 全局变量变化处 */ 794 AddPosPart((x1 + x2 + x2) / 3, (y1 + y1 + y2) / 3, (x2 - x1)*(y2 - y1) / 2, xtr, ytr, wtr); 795 // triangle 全局变量变化处 796 } 797 else 798 { 799 AddNegPart((x2 + x1) / 2, y1 / 2, (x2 - x1) * y1, xtl, ytl, wtl); 800 // rectangle 全局变量变化处 801 AddNegPart((x1 + x2 + x2) / 3, (y1 + y1 + y2) / 3, (x2 - x1)*(y2 - y1) / 2, xtl, ytl, wtl); 802 // triangle 全局变量变化处 803 } 804 } 805 POINT cg_simple(int vcount, POINT polygon[]) 806 { 807 double xtr, ytr, wtr, xtl, ytl, wtl; 808 //注意： 如果把xtr,ytr,wtr,xtl,ytl,wtl改成全局变量后这里要删去 809 POINT p1, p2, tp; 810 xtr = ytr = wtr = 0.0; 811 xtl = ytl = wtl = 0.0; 812 for (int i = 0; i<vcount; i++) 813 { 814 p1 = polygon[i]; 815 p2 = polygon[(i + 1) % vcount]; 816 AddRegion(p1.x, p1.y, p2.x, p2.y, xtr, ytr, wtr, xtl, ytl, wtl); //全局变量变化处 817 } 818 tp.x = (wtr*xtr + wtl*xtl) / (wtr + wtl); 819 tp.y = (wtr*ytr + wtl*ytl) / (wtr + wtl); 820 return tp; 821 } 822 // 求凸多边形的重心,要求输入多边形按逆时针排序 823 POINT gravitycenter(int vcount, POINT polygon[]) 824 { 825 POINT tp; 826 double x, y, s, x0, y0, cs, k; 827 x = 0; y = 0; s = 0; 828 for (int i = 1; i<vcount - 1; i++) 829 { 830 x0 = (polygon[0].x + polygon[i].x + polygon[i + 1].x) / 3; 831 y0 = (polygon[0].y + polygon[i].y + polygon[i + 1].y) / 3; //求当前三角形的重心 832 cs = multiply(polygon[i], polygon[i + 1], polygon[0]) / 2; 833 //三角形面积可以直接利用该公式求解 834 if (abs(s)<1e-20) 835 { 836 x = x0; y = y0; s += cs; continue; 837 } 838 k = cs / s; //求面积比例 839 x = (x + k*x0) / (1 + k); 840 y = (y + k*y0) / (1 + k); 841 s += cs; 842 } 843 tp.x = x; 844 tp.y = y; 845 return tp; 846 } 847 848 /************************************************ 849 给定一简单多边形，找出一个肯定在该多边形内的点 850 定理1 ：每个多边形至少有一个凸顶点 851 定理2 ：顶点数>=4的简单多边形至少有一条对角线 852 结论 ： x坐标最大，最小的点肯定是凸顶点 853 y坐标最大，最小的点肯定是凸顶点 854 ************************************************/ 855 POINT a_point_insidepoly(int vcount, POINT polygon[]) 856 { 857 POINT v, a, b, r; 858 int i, index; 859 v = polygon[0]; 860 index = 0; 861 for (i = 1; i<vcount; i++) //寻找一个凸顶点 862 { 863 if (polygon[i].y<v.y) 864 { 865 v = polygon[i]; 866 index = i; 867 } 868 } 869 a = polygon[(index - 1 + vcount) % vcount]; //得到v的前一个顶点 870 b = polygon[(index + 1) % vcount]; //得到v的后一个顶点 871 POINT tri[3], q; 872 tri[0] = a; tri[1] = v; tri[2] = b; 873 double md = INF; 874 int in1 = index; 875 bool bin = false; 876 for (i = 0; i<vcount; i++) //寻找在三角形avb内且离顶点v最近的顶点q 877 { 878 if (i == index)continue; 879 if (i == (index - 1 + vcount) % vcount)continue; 880 if (i == (index + 1) % vcount)continue; 881 if (!InsideConvexPolygon(3, tri, polygon[i]))continue; 882 bin = true; 883 if (dist(v, polygon[i])<md) 884 { 885 q = polygon[i]; 886 md = dist(v, q); 887 } 888 } 889 if (!bin) //没有顶点在三角形avb内，返回线段ab中点 890 { 891 r.x = (a.x + b.x) / 2; 892 r.y = (a.y + b.y) / 2; 893 return r; 894 } 895 r.x = (v.x + q.x) / 2; //返回线段vq的中点 896 r.y = (v.y + q.y) / 2; 897 return r; 898 } 899 /*********************************************************************************************** 900 求从多边形外一点p出发到一个简单多边形的切线,如果存在返回切点,其中rp点是右切点,lp是左切点 901 注意：p点一定要在多边形外 ,输入顶点序列是逆时针排列 902 原 理： 如果点在多边形内肯定无切线;凸多边形有唯一的两个切点,凹多边形就可能有多于两个的切点; 903 如果polygon是凸多边形，切点只有两个只要找到就可以,可以化简此算法 904 如果是凹多边形还有一种算法可以求解:先求凹多边形的凸包,然后求凸包的切线 905 /***********************************************************************************************/ 906 void pointtangentpoly(int vcount, POINT polygon[], POINT p, POINT &rp, POINT &lp) 907 { 908 LINESEG ep, en; 909 bool blp, bln; 910 rp = polygon[0]; 911 lp = polygon[0]; 912 for (int i = 1; i<vcount; i++) 913 { 914 ep.s = polygon[(i + vcount - 1) % vcount]; 915 ep.e = polygon[i]; 916 en.s = polygon[i]; 917 en.e = polygon[(i + 1) % vcount]; 918 blp = multiply(ep.e, p, ep.s) >= 0; // p is to the left of pre edge 919 bln = multiply(en.e, p, en.s) >= 0; // p is to the left of next edge 920 if (!blp&&bln) 921 { 922 if (multiply(polygon[i], rp, p)>0) // polygon[i] is above rp 923 rp = polygon[i]; 924 } 925 if (blp && !bln) 926 { 927 if (multiply(lp, polygon[i], p)>0) // polygon[i] is below lp 928 lp = polygon[i]; 929 } 930 } 931 return; 932 } 933 // 如果多边形polygon的核存在，返回true，返回核上的一点p.顶点按逆时针方向输入 934 bool core_exist(int vcount, POINT polygon[], POINT &p) 935 { 936 int i, j, k; 937 LINESEG l; 938 LINE lineset[MAXV]; 939 for (i = 0; i<vcount; i++) 940 { 941 lineset[i] = makeline(polygon[i], polygon[(i + 1) % vcount]); 942 } 943 for (i = 0; i<vcount; i++) 944 { 945 for (j = 0; j<vcount; j++) 946 { 947 if (i == j)continue; 948 if (lineintersect(lineset[i], lineset[j], p)) 949 { 950 for (k = 0; k<vcount; k++) 951 { 952 l.s = polygon[k]; 953 l.e = polygon[(k + 1) % vcount]; 954 if (multiply(p, l.e, l.s)>0) 955 //多边形顶点按逆时针方向排列，核肯定在每条边的左侧或边上 956 break; 957 } 958 if (k == vcount) //找到了一个核上的点 959 break; 960 } 961 } 962 if (j<vcount) break; 963 } 964 if (i<vcount) 965 return true; 966 else 967 return false; 968 } 969 /************************* 970 * * 971 * 圆的基本运算 * 972 * * 973 \*************************/ 974 /****************************************************************************** 975 返回值 ： 点p在圆内(包括边界)时，返回true 976 用途 ： 因为圆为凸集，所以判断点集，折线，多边形是否在圆内时， 977 只需要逐一判断点是否在圆内即可。 978 *******************************************************************************/ 979 bool point_in_circle(POINT o, double r, POINT p) 980 { 981 double d2 = (p.x - o.x)*(p.x - o.x) + (p.y - o.y)*(p.y - o.y); 982 double r2 = r*r; 983 return d2<r2 || abs(d2 - r2)<EP; 984 } 985 /****************************************************************************** 986 用 途 ：求不共线的三点确定一个圆 987 输 入 ：三个点p1,p2,p3 988 返回值 ：如果三点共线，返回false；反之，返回true。圆心由q返回，半径由r返回 989 *******************************************************************************/ 990 bool cocircle(POINT p1, POINT p2, POINT p3, POINT &q, double &r) 991 { 992 double x12 = p2.x - p1.x; 993 double y12 = p2.y - p1.y; 994 double x13 = p3.x - p1.x; 995 double y13 = p3.y - p1.y; 996 double z2 = x12*(p1.x + p2.x) + y12*(p1.y + p2.y); 997 double z3 = x13*(p1.x + p3.x) + y13*(p1.y + p3.y); 998 double d = 2.0*(x12*(p3.y - p2.y) - y12*(p3.x - p2.x)); 999 if (abs(d)<EP) //共线，圆不存在 1000 return false; 1001 q.x = (y13*z2 - y12*z3) / d; 1002 q.y = (x12*z3 - x13*z2) / d; 1003 r = dist(p1, q); 1004 return true; 1005 } 1006 int line_circle(LINE l, POINT o, double r, POINT &p1, POINT &p2) 1007 { 1008 return true; 1009 } 1010 1011 /**************************1012 * * 1013 * 矩形的基本运算 * 1014 * * 1015 \**************************/ 1016 /* 1017 说明：因为矩形的特殊性，常用算法可以化简： 1018 1.判断矩形是否包含点 1019 只要判断该点的横坐标和纵坐标是否夹在矩形的左右边和上下边之间。 1020 2.判断线段、折线、多边形是否在矩形中 1021 因为矩形是个凸集，所以只要判断所有端点是否都在矩形中就可以了。 1022 3.判断圆是否在矩形中 1023 圆在矩形中的充要条件是：圆心在矩形中且圆的半径小于等于圆心到矩形四边的距离的最小值。 1024 */ 1025 // 已知矩形的三个顶点(a,b,c)，计算第四个顶点d的坐标. 注意：已知的三个顶点可以是无序的 1026 POINT rect4th(POINT a, POINT b, POINT c) 1027 { 1028 POINT d; 1029 if (abs(dotmultiply(a, b, c))<EP) // 说明c点是直角拐角处 1030 { 1031 d.x = a.x + b.x - c.x; 1032 d.y = a.y + b.y - c.y; 1033 } 1034 if (abs(dotmultiply(a, c, b))<EP) // 说明b点是直角拐角处 1035 { 1036 d.x = a.x + c.x - b.x; 1037 d.y = a.y + c.y - b.x; 1038 } 1039 if (abs(dotmultiply(c, b, a))<EP) // 说明a点是直角拐角处 1040 { 1041 d.x = c.x + b.x - a.x; 1042 d.y = c.y + b.y - a.y; 1043 } 1044 return d; 1045 } 1046 1047 /*************************1048 * * 1049 * 常用算法的描述 * 1050 * * 1051 \*************************/ 1052 /* 1053 尚未实现的算法： 1054 1. 求包含点集的最小圆 1055 2. 求多边形的交 1056 3. 简单多边形的三角剖分 1057 4. 寻找包含点集的最小矩形 1058 5. 折线的化简 1059 6. 判断矩形是否在矩形中 1060 7. 判断矩形能否放在矩形中 1061 8. 矩形并的面积与周长 1062 9. 矩形并的轮廓 1063 10.矩形并的闭包 1064 11.矩形的交 1065 12.点集中的最近点对 1066 13.多边形的并 1067 14.圆的交与并 1068 15.直线与圆的关系 1069 16.线段与圆的关系 1070 17.求多边形的核监视摄象机 1071 18.求点集中不相交点对 railwai 1072 *//* 1073 寻找包含点集的最小矩形 1074 原理：该矩形至少一条边与点集的凸壳的某条边共线 1075 First take the convex hull of the points. Let the resulting convex 1076 polygon be P. It has been known for some time that the minimum 1077 area rectangle enclosing P must have one rectangle side flush with 1078 (i.e., collinear with and overlapping) one edge of P. This geometric 1079 fact was used by Godfried Toussaint to develop the "rotating calipers" 1080 algorithm in a hard-to-find 1983 paper, "Solving Geometric Problems 1081 with the Rotating Calipers" (Proc. IEEE MELECON). The algorithm 1082 rotates a surrounding rectangle from one flush edge to the next, 1083 keeping track of the minimum area for each edge. It achieves O(n) 1084 time (after hull computation). See the "Rotating Calipers Homepage" 1085 http://www.cs.mcgill.ca/~orm/rotcal.frame.html for a description 1086 and applet. 1087 *//* 1088 折线的化简 伪码如下： 1089 Input: tol = the approximation tolerance 1090 L = {V0,V1,,Vn-1} is any n-vertex polyline 1091 1092 Set start = 0; 1093 Set k = 0; 1094 Set W0 = V0; 1095 for each vertex Vi (i=1,n-1) 1096 { 1097 if Vi is within tol from Vstart 1098 then ignore it, and continue with the next vertex 1099 1100 Vi is further than tol away from Vstart 1101 so add it as a new vertex of the reduced polyline 1102 Increment k++; 1103 Set Wk = Vi; 1104 Set start = i; as the new initial vertex 1105 } 1106 1107 Output: W = {W0,W1,,Wk-1} = the k-vertex simplified polyline 1108 */ 1109 /********************1110 * * 1111 * 补充 * 1112 * * 1113 \********************/ 1114 1115 //两圆关系： 1116 /* 两圆： 1117 相离： return 1； 1118 外切： return 2； 1119 相交： return 3； 1120 内切： return 4； 1121 内含： return 5； 1122 */ 1123 int CircleRelation(POINT p1, double r1, POINT p2, double r2) 1124 { 1125 double d = sqrt((p1.x - p2.x)*(p1.x - p2.x) + (p1.y - p2.y)*(p1.y - p2.y)); 1126 1127 if (fabs(d - r1 - r2) < EP) // 必须保证前两个if先被判定！ 1128 return 2; 1129 if (fabs(d - fabs(r1 - r2)) < EP) 1130 return 4; 1131 if (d > r1 + r2) 1132 return 1; 1133 if (d < fabs(r1 - r2)) 1134 return 5; 1135 if (fabs(r1 - r2) < d && d < r1 + r2) 1136 return 3; 1137 return 0; // indicate an error! 1138 } 1139 //判断圆是否在矩形内： 1140 // 判定圆是否在矩形内，是就返回true（设矩形水平，且其四个顶点由左上开始按顺时针排列） 1141 // 调用ptoldist函数，在第4页 1142 bool CircleRecRelation(POINT pc, double r, POINT pr1, POINT pr2, POINT pr3, POINT pr4) 1143 { 1144 if (pr1.x < pc.x && pc.x < pr2.x && pr3.y < pc.y && pc.y < pr2.y) 1145 { 1146 LINESEG line1(pr1, pr2); 1147 LINESEG line2(pr2, pr3); 1148 LINESEG line3(pr3, pr4); 1149 LINESEG line4(pr4, pr1); 1150 if (r<ptoldist(pc, line1) && r<ptoldist(pc, line2) && r<ptoldist(pc, line3) && r<ptoldist(pc, line4)) 1151 return true; 1152 } 1153 return false; 1154 } 1155 //点到平面的距离： 1156 //点到平面的距离,平面用一般式表示ax+by+cz+d=0 1157 double P2planeDist(double x, double y, double z, double a, double b, double c, double d) 1158 { 1159 return fabs(a*x + b*y + c*z + d) / sqrt(a*a + b*b + c*c); 1160 } 1161 //点是否在直线同侧： 1162 //两个点是否在直线同侧，是则返回true 1163 bool SameSide(POINT p1, POINT p2, LINE line) 1164 { 1165 return (line.a * p1.x + line.b * p1.y + line.c) * 1166 (line.a * p2.x + line.b * p2.y + line.c) > 0; 1167 } 1168 //镜面反射线： 1169 // 已知入射线、镜面，求反射线。 1170 // a1,b1,c1为镜面直线方程(a1 x + b1 y + c1 = 0 ,下同)系数; 1171 //a2,b2,c2为入射光直线方程系数; 1172 //a,b,c为反射光直线方程系数. 1173 // 光是有方向的，使用时注意：入射光向量:<-b2,a2>；反射光向量:<b,-a>. 1174 // 不要忘记结果中可能会有"negative zeros" 1175 void reflect(double a1, double b1, double c1, double a2, double b2, double c2, double &a, double &b, double &c) 1176 { 1177 double n, m; 1178 double tpb, tpa; 1179 tpb = b1*b2 + a1*a2; 1180 tpa = a2*b1 - a1*b2; 1181 m = (tpb*b1 + tpa*a1) / (b1*b1 + a1*a1); 1182 n = (tpa*b1 - tpb*a1) / (b1*b1 + a1*a1); 1183 if (fabs(a1*b2 - a2*b1)<1e-20) 1184 { 1185 a = a2; b = b2; c = c2; 1186 return; 1187 } 1188 double xx, yy; //(xx,yy)是入射线与镜面的交点。 1189 xx = (b1*c2 - b2*c1) / (a1*b2 - a2*b1); 1190 yy = (a2*c1 - a1*c2) / (a1*b2 - a2*b1); 1191 a = n; 1192 b = -m; 1193 c = m*yy - xx*n; 1194 } 1195 //矩形包含： 1196 // 矩形2（C，D）是否在1（A，B）内 1197 bool r2inr1(double A, double B, double C, double D) 1198 { 1199 double X, Y, L, K, DMax; 1200 if (A < B) 1201 { 1202 double tmp = A; 1203 A = B; 1204 B = tmp; 1205 } 1206 if (C < D) 1207 { 1208 double tmp = C; 1209 C = D; 1210 D = tmp; 1211 } 1212 if (A > C && B > D) // trivial case 1213 return true; 1214 else 1215 if (D >= B) 1216 return false; 1217 else 1218 { 1219 X = sqrt(A * A + B * B); // outer rectangle‘s diagonal 1220 Y = sqrt(C * C + D * D); // inner rectangle‘s diagonal 1221 if (Y < B) // check for marginal conditions 1222 return true; // the inner rectangle can freely rotate inside 1223 else 1224 if (Y > X) 1225 return false; 1226 else 1227 { 1228 L = (B - sqrt(Y * Y - A * A)) / 2; 1229 K = (A - sqrt(Y * Y - B * B)) / 2; 1230 DMax = sqrt(L * L + K * K); 1231 if (D >= DMax) 1232 return false; 1233 else 1234 return true; 1235 } 1236 } 1237 } 1238 //两圆交点： 1239 // 两圆已经相交（相切） 1240 void c2point(POINT p1, double r1, POINT p2, double r2, POINT &rp1, POINT &rp2) 1241 { 1242 double a, b, r; 1243 a = p2.x - p1.x; 1244 b = p2.y - p1.y; 1245 r = (a*a + b*b + r1*r1 - r2*r2) / 2; 1246 if (a == 0 && b != 0) 1247 { 1248 rp1.y = rp2.y = r / b; 1249 rp1.x = sqrt(r1*r1 - rp1.y*rp1.y); 1250 rp2.x = -rp1.x; 1251 } 1252 else if (a != 0 && b == 0) 1253 { 1254 rp1.x = rp2.x = r / a; 1255 rp1.y = sqrt(r1*r1 - rp1.x*rp2.x); 1256 rp2.y = -rp1.y; 1257 } 1258 else if (a != 0 && b != 0) 1259 { 1260 double delta; 1261 delta = b*b*r*r - (a*a + b*b)*(r*r - r1*r1*a*a); 1262 rp1.y = (b*r + sqrt(delta)) / (a*a + b*b); 1263 rp2.y = (b*r - sqrt(delta)) / (a*a + b*b); 1264 rp1.x = (r - b*rp1.y) / a; 1265 rp2.x = (r - b*rp2.y) / a; 1266 } 1267 1268 rp1.x += p1.x; 1269 rp1.y += p1.y; 1270 rp2.x += p1.x; 1271 rp2.y += p1.y; 1272 } 1273 //两圆公共面积： 1274 // 必须保证相交 1275 double c2area(POINT p1, double r1, POINT p2, double r2) 1276 { 1277 POINT rp1, rp2; 1278 c2point(p1, r1, p2, r2, rp1, rp2); 1279 1280 if (r1>r2) //保证r2>r1 1281 { 1282 swap(p1, p2); 1283 swap(r1, r2); 1284 } 1285 double a, b, rr; 1286 a = p1.x - p2.x; 1287 b = p1.y - p2.y; 1288 rr = sqrt(a*a + b*b); 1289 1290 double dx1, dy1, dx2, dy2; 1291 double sita1, sita2; 1292 dx1 = rp1.x - p1.x; 1293 dy1 = rp1.y - p1.y; 1294 dx2 = rp2.x - p1.x; 1295 dy2 = rp2.y - p1.y; 1296 sita1 = acos((dx1*dx2 + dy1*dy2) / r1 / r1); 1297 1298 dx1 = rp1.x - p2.x; 1299 dy1 = rp1.y - p2.y; 1300 dx2 = rp2.x - p2.x; 1301 dy2 = rp2.y - p2.y; 1302 sita2 = acos((dx1*dx2 + dy1*dy2) / r2 / r2); 1303 double s = 0; 1304 if (rr<r2)//相交弧为优弧 1305 s = r1*r1*(PI - sita1 / 2 + sin(sita1) / 2) + r2*r2*(sita2 - sin(sita2)) / 2; 1306 else//相交弧为劣弧 1307 s = (r1*r1*(sita1 - sin(sita1)) + r2*r2*(sita2 - sin(sita2))) / 2; 1308 1309 return s; 1310 } 1311 //圆和直线关系： 1312 //0----相离 1----相切 2----相交 1313 int clpoint(POINT p, double r, double a, double b, double c, POINT &rp1, POINT &rp2) 1314 { 1315 int res = 0; 1316 1317 c = c + a*p.x + b*p.y; 1318 double tmp; 1319 if (a == 0 && b != 0) 1320 { 1321 tmp = -c / b; 1322 if (r*r<tmp*tmp) 1323 res = 0; 1324 else if (r*r == tmp*tmp) 1325 { 1326 res = 1; 1327 rp1.y = tmp; 1328 rp1.x = 0; 1329 } 1330 else 1331 { 1332 res = 2; 1333 rp1.y = rp2.y = tmp; 1334 rp1.x = sqrt(r*r - tmp*tmp); 1335 rp2.x = -rp1.x; 1336 } 1337 } 1338 else if (a != 0 && b == 0) 1339 { 1340 tmp = -c / a; 1341 if (r*r<tmp*tmp) 1342 res = 0; 1343 else if (r*r == tmp*tmp) 1344 { 1345 res = 1; 1346 rp1.x = tmp; 1347 rp1.y = 0; 1348 } 1349 else 1350 { 1351 res = 2; 1352 rp1.x = rp2.x = tmp; 1353 rp1.y = sqrt(r*r - tmp*tmp); 1354 rp2.y = -rp1.y; 1355 } 1356 } 1357 else if (a != 0 && b != 0) 1358 { 1359 double delta; 1360 delta = b*b*c*c - (a*a + b*b)*(c*c - a*a*r*r); 1361 if (delta<0) 1362 res = 0; 1363 else if (delta == 0) 1364 { 1365 res = 1; 1366 rp1.y = -b*c / (a*a + b*b); 1367 rp1.x = (-c - b*rp1.y) / a; 1368 } 1369 else 1370 { 1371 res = 2; 1372 rp1.y = (-b*c + sqrt(delta)) / (a*a + b*b); 1373 rp2.y = (-b*c - sqrt(delta)) / (a*a + b*b); 1374 rp1.x = (-c - b*rp1.y) / a; 1375 rp2.x = (-c - b*rp2.y) / a; 1376 } 1377 } 1378 rp1.x += p.x; 1379 rp1.y += p.y; 1380 rp2.x += p.x; 1381 rp2.y += p.y; 1382 return res; 1383 } 1384 //内切圆： 1385 void incircle(POINT p1, POINT p2, POINT p3, POINT &rp, double &r) 1386 { 1387 double dx31, dy31, dx21, dy21, d31, d21, a1, b1, c1; 1388 dx31 = p3.x - p1.x; 1389 dy31 = p3.y - p1.y; 1390 dx21 = p2.x - p1.x; 1391 dy21 = p2.y - p1.y; 1392 1393 d31 = sqrt(dx31*dx31 + dy31*dy31); 1394 d21 = sqrt(dx21*dx21 + dy21*dy21); 1395 a1 = dx31*d21 - dx21*d31; 1396 b1 = dy31*d21 - dy21*d31; 1397 c1 = a1*p1.x + b1*p1.y; 1398 1399 double dx32, dy32, dx12, dy12, d32, d12, a2, b2, c2; 1400 dx32 = p3.x - p2.x; 1401 dy32 = p3.y - p2.y; 1402 dx12 = -dx21; 1403 dy12 = -dy21; 1404 1405 d32 = sqrt(dx32*dx32 + dy32*dy32); 1406 d12 = d21; 1407 a2 = dx12*d32 - dx32*d12; 1408 b2 = dy12*d32 - dy32*d12; 1409 c2 = a2*p2.x + b2*p2.y; 1410 1411 rp.x = (c1*b2 - c2*b1) / (a1*b2 - a2*b1); 1412 rp.y = (c2*a1 - c1*a2) / (a1*b2 - a2*b1); 1413 r = fabs(dy21*rp.x - dx21*rp.y + dx21*p1.y - dy21*p1.x) / d21; 1414 } 1415 //求切点： 1416 // p---圆心坐标， r---圆半径， sp---圆外一点， rp1,rp2---切点坐标 1417 void cutpoint(POINT p, double r, POINT sp, POINT &rp1, POINT &rp2) 1418 { 1419 POINT p2; 1420 p2.x = (p.x + sp.x) / 2; 1421 p2.y = (p.y + sp.y) / 2; 1422 1423 double dx2, dy2, r2; 1424 dx2 = p2.x - p.x; 1425 dy2 = p2.y - p.y; 1426 r2 = sqrt(dx2*dx2 + dy2*dy2); 1427 c2point(p, r, p2, r2, rp1, rp2); 1428 } 1429 //线段的左右旋： 1430 /* l2在l1的左/右方向（l1为基准线） 1431 返回 0 ： 重合； 1432 返回 1 ： 右旋； 1433 返回 –1 ： 左旋； 1434 */ 1435 int rotat(LINESEG l1, LINESEG l2) 1436 { 1437 double dx1, dx2, dy1, dy2; 1438 dx1 = l1.s.x - l1.e.x; 1439 dy1 = l1.s.y - l1.e.y; 1440 dx2 = l2.s.x - l2.e.x; 1441 dy2 = l2.s.y - l2.e.y; 1442 1443 double d; 1444 d = dx1*dy2 - dx2*dy1; 1445 if (d == 0) 1446 return 0; 1447 else if (d>0) 1448 return -1; 1449 else 1450 return 1; 1451 } 1452 /* 1453 公式： 1454 1455 球坐标公式： 1456 直角坐标为 P(x, y, z) 时，对应的球坐标是(rsinφcosθ, rsinφsinθ, rcosφ),其中φ是向量OP与Z轴的夹角，范围[0，π]；是OP在XOY面上的投影到X轴的旋角，范围[0，2π] 1457 1458 直线的一般方程转化成向量方程： 1459 ax+by+c=0 1460 x-x0 y-y0 1461 ------ = ------- // (x0,y0)为直线上一点，m,n为向量 1462 m n 1463 转换关系： 1464 a=n；b=-m；c=m·y0-n·x0； 1465 m=-b; n=a; 1466 1467 三点平面方程： 1468 三点为P1，P2，P3 1469 设向量 M1=P2-P1; M2=P3-P1; 1470 平面法向量： M=M1 x M2 （） 1471 平面方程： M.i(x-P1.x)+M.j(y-P1.y)+M.k(z-P1.z)=0 1472 */
temp、
声明：以上内容来自用户投稿及互联网公开渠道收集整理发布，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任，若内容有误或涉及侵权可进行投诉： 投诉/举报 工作人员会在5个工作日内联系你，一经查实，本站将立刻删除涉嫌侵权内容。