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