简单的计算几何
已解决问题：
判断点是否在线段上
判断两线段是否相交
判断点是否在多边形内
判断线段、折线、多边形是否在多变形内
判断上述是否在圆内
计算上述与线段及直线的交点
凸包
待解决：
半平面交
1 #include<ctime> 2 #include<cstdio> 3 #include<cstring> 4 #include<cstdlib> 5 #include<algorithm> 6 #include<cmath> 7 #include<vector> 8 #include<stack> 9 #include<queue> 10 11 using namespace std; 12 //================================ struct ================================= 13 struct V{ 14 double px,py; 15 V operator - (const V c){ 16 V a; 17 a.px = px - c.px,a.py = py - c.py; 18 return a; 19 } 20 V operator + (const V c){ 21 V a; 22 a.px = px + c.px,a.py = py + c.py; 23 return a; 24 } 25 V operator * (const double c){ 26 V a = *this; 27 a.px *= c,a.py *= c; 28 return a; 29 } 30 V operator / (const double c){ 31 V a = *this; 32 a.px /= c,a.py /= c; 33 return a; 34 } 35 bool operator <= (const V c){ 36 //return px <= c.px && py <= c.py; 37 return px*px + py*py <= c.px*c.px + c.py*c.py; 38 } 39 bool operator == (const V c){ 40 return px == c.px && py == c.py; 41 } 42 double MO2(){ 43 return px * px + py * py; 44 } 45 double DIS(V p1,V p2,int pa){ 46 if(pa == 1){ 47 return abs(py - p1.py); 48 } 49 if(pa == 2){ 50 return abs(px - p1.px); 51 } 52 if(pa == 0){ 53 double k = (p2.py-p1.py)/(p2.px-p1.px),x0,y0; 54 x0 = (k*k * p1.px + k * (py - p1.py) + px) / (k*k + 1); 55 y0 = k * (px - p1.px) + p1.py; 56 return sqrt((px-x0)*(px-x0) + (py-y0)*(py-y0)); 57 } 58 } 59 double CP(V c){ 60 return px*c.py - py*c.px; 61 }//cross product 62 double DP(V c){ 63 return px*c.px + py*c.py; 64 }//dot product 65 bool min(V a,V b){ 66 if(b <= a) swap(a,b); 67 return a <= *this; 68 } 69 bool max(V a,V b){ 70 if(a <= b) swap(a,b); 71 return *this <= a; 72 } 73 };//vector 74 75 struct L{ 76 V p1,p2; 77 int p0; 78 /* 79 switch(p0){ 80 case 1: beeline 81 case 2: segment 82 case 3: ray 83 } 84 */ 85 int pa; 86 /* 87 switch(pa){ 88 case 1: parallel-y; 89 case 2: parallel-x; 90 case 0: k < oo && k != 0; 91 } 92 */ 93 double k; 94 void getk(){ 95 if(p1.px == p2.px) pa = 1; 96 else if(p1.py == p2.py) pa = 2; 97 else pa = 0,k = (p1.py-p2.py) / (p1.px-p2.px); 98 } 99 bool IP(V q){ 100 return q.min(p1,p2) && q.max(p1,p2); 101 }//in the ploygon of the line 102 bool PA(L c){ 103 if(pa != c.pa) return false; 104 if(pa) return true; 105 else if(k == c.k) return true; 106 else return false; 107 } 108 bool CR(L c){ 109 if((*this).PA(c)) return false; 110 V v1 = p2 - p1,v2 = c.p1 - p1,v3 = c.p2 - p1; 111 double exp = v1.CP(v2) * v1.CP(v3); 112 //if(v1.CP(v2) * v1.CP(v3) < 0) return true; 113 switch(p0){ 114 case 1://B_S 115 return exp <= 0; 116 break; 117 case 2://S_S 118 return exp < 0; 119 break; 120 case 3://R_S 121 if(exp > 0) return false; 122 else if(!exp){ 123 V v = v1.CP(v2)? c.p2:c.p1; 124 if(p1-v <= p2-v) return false; 125 else return true; 126 } 127 else{ 128 V v4 = c.p1 - p2,v5 = c.p2 - p2; 129 double s1 = abs(v2.CP(v3)),s2 = abs(v4.CP(v5)); 130 return s1 > s2; 131 } 132 break; 133 } 134 }//line cross 135 };//line 136 137 struct BL{ 138 vector<L> E; 139 };//broken line 140 141 struct PL{ 142 vector<L> E; 143 };//ploygon 144 145 struct O{ 146 V o; 147 double r; 148 };//circle 149 150 //================================ data =================================== 151 double CP(V a,V b){ 152 return a.px*b.py - a.py*b.px; 153 }//cross product 154 155 double DP(V a,V b){ 156 return a.px*b.px + a.py*b.py; 157 }//dot product 158 159 double P_P_DIST(V a,V b){ 160 return sqrt((a.px-b.px)*(a.px-b.px) + (a.py-b.py)*(a.py-b.py)); 161 }//disdance for points 162 163 double P_L_DIST(V a,L b){ 164 double ret = a.DIS(b.p1,b.p2,b.pa); 165 if(b.p0 != 1){ 166 ret = min(ret,P_P_DIST(a,b.p1)); 167 } 168 if(b.p0 == 2){ 169 ret = min(ret,P_P_DIST(a,b.p2)); 170 } 171 return ret; 172 }//least distance for point to line(segment/beeline/ray) 173 174 double P_BL_DIST(V p,BL bl){ 175 double ret = oo; 176 for(int i = 0;i < bl.E.size();++ i){ 177 ret = min(ret,P_L_DIST(p,bl.E[i])); 178 } 179 return ret; 180 }//least distance for point to broken line 181 182 double P_PL_DIST(V p,PL pl){ 183 double ret = oo; 184 for(int i = 0;i < pl.E.size();++ i){ 185 ret = min(ret,P_L_DIST(p,pl.E[i])); 186 } 187 return ret; 188 }//least distance for point to ploygon 189 190 double P_O_DIST(V p,O o){ 191 return o.r - P_P_DIST(p,o.o); 192 }//least distance for point to circle 193 194 double AREA(PL pl){ 195 V v0 = pl.E[0].p1; 196 double ret = 0; 197 for(int i = 0;i < pl.E.size();++ i){ 198 ret += v0.CP(pl.E[i].p2-pl.E[i].p1); 199 v0 = v0 + (pl.E[i].p2-pl.E[i].p1); 200 } 201 return ret; 202 }//the measure of area 203 204 //================================= judge ================================= 205 bool P_L_ON(V q,L p){ 206 /*double x1 = p.p1.px,x2 = p.p2.px,y1 = p.p1.py,y2 = p.p2.py, 207 x0 = q.px,y0 = q.py; 208 if((!p.p0) && (!(x0 <= max(x1,x2) && x0 >= min (x1,x2) && y0 <= max(y1,y2) && y0 >= min(y1,y2)))) 209 return false;*/ 210 //if((!p.p0) && (!(q.min(p.p1,p.p2) && q.max(p.p1,p.p2)))) return false; 211 if((p.p0 == 2) && (!p.IP(q))) return false; 212 int pd = CP(p.p1-q,p.p2-q); 213 if(pd != 0) return false; 214 pd = DP(p.p1-q,p.p2-q); 215 if(pd <= 0) return true; 216 else{ 217 if(p.p0 == 1) return true; 218 if(p.p0 == 2) return false; 219 if(p.p0 == 3){ 220 if(p.p1-q <= p.p2-q) return false; 221 else return true; 222 } 223 } 224 //if(!CP(q-p.p1,p.p2-p.p1)) return true; 225 //else return false; 226 }//wheather a point is on a line(beeline/segment/ray) 227 228 bool S_S_CR(L a,L b){ 229 return (a.CR(b) && b.CR(a)) || ((P_L_ON(a.p1,b) || P_L_ON(a.p2,b) || P_L_ON(b.p1,a) || P_L_ON(b.p2,a))/* && ((a.p2-a.p1).DP(b.p1-a.p1) * (a.p2-a.p1).DP(b.p2-a.p1) == 0)*/); 230 }//segments' cross 231 232 bool S_L_CR(L a,L b){ 233 return b.CR(a); 234 }//segment and beeline's cross 235 236 bool S_R_CR(L a,L b){ 237 return b.CR(a); 238 }//segment and ray's cross 239 240 bool P_PL_IN(V p,PL pl){ 241 L o; 242 o.p1 = p,o.p2.px = p.px-1,o.p2.py = p.py,o.p0 = 3; 243 int cnt = 0; 244 for(int i = 0;i < pl.E.size();++ i){ 245 if(P_L_ON(p,pl.E[i])) return true; 246 if(pl.E[i].pa != 1){ 247 if(P_L_ON(pl.E[i].p1,o) || (P_L_ON(pl.E[i].p2,o))){ 248 V v1 = P_L_ON(pl.E[i].p1,o)? pl.E[i].p1:pl.E[i].p2, 249 v2 = P_L_ON(pl.E[i].p1,o)? pl.E[i].p2:pl.E[i].p1; 250 if(v1.py > v2.py) cnt ++; 251 } 252 } 253 else if(S_R_CR(pl.E[i],o)) cnt ++; 254 } 255 return cnt & 1; 256 }//point in ploygon 257 258 bool S_PL_IN(L p,PL pl){ 259 if(!(P_PL_IN(p.p1,pl) && P_PL_IN(p.p2,pl))) return false; 260 vector<V> v; 261 for(int i = 0;i < pl.E.size();++ i){ 262 if(P_L_ON(p.p1,pl.E[i]) || P_L_ON(p.p2,pl.E[i])){ 263 if(P_L_ON(p.p1,pl.E[i])) v.push_back(p.p1); 264 if(P_L_ON(p.p2,pl.E[i])) v.push_back(p.p2); 265 } 266 else if(P_L_ON(pl.E[i].p1,p) || P_L_ON(pl.E[i].p2,p)){ 267 if(P_L_ON(pl.E[i].p1,p)) v.push_back(pl.E[i].p1); 268 if(P_L_ON(pl.E[i].p2,p)) v.push_back(pl.E[i].p2); 269 } 270 else if(p.CR(pl.E[i])) return false; 271 } 272 for(int i = 0;i < v.size() - 1;++ i){ 273 if(!P_PL_IN((v[i] + v[i+1])/2,pl)) return false; 274 } 275 return true; 276 }//segment in ploygon 277 278 bool BL_PL_IN(BL bl,PL pl){ 279 for(int i = 0;i < bl.E.size();++ i){ 280 if(!S_PL_IN(bl.E[i],pl)) return false; 281 } 282 return true; 283 }//broken line in ploygon 284 285 bool PL_PL_IN(PL pl1,PL pl2){ 286 for(int i = 0;i < pl1.E.size();++ i){ 287 if(!S_PL_IN(pl1.E[i],pl2)) return false; 288 } 289 return true; 290 }//ploygon in ploygon 291 292 bool O_PL_IN(O o,PL pl){ 293 double dis = oo; 294 for(int i = 0;i < pl.E.size();++ i){ 295 dis = min(dis,P_L_DIST(o.o,pl.E[i])); 296 } 297 return o.r <= dis; 298 }//circle in ploygon 299 300 bool P_O_IN(V p,O o){ 301 return (p-o.o).MO2() <= o.r * o.r; 302 }//point in circle 303 304 bool S_O_IN(L s,O o){ 305 return P_O_IN(s.p1,o) && P_O_IN(s.p2,o); 306 }//segment in circle 307 308 bool BL_O_IN(BL bl,O o){ 309 for(int i = 0;i < bl.E.size();++ i){ 310 if(!P_O_IN(bl.E[i].p1,o)) return false; 311 } 312 return P_O_IN(bl.E[bl.E.size()-1].p2,o); 313 }//broken line in circle 314 315 bool PL_O_IN(PL pl,O o){ 316 for(int i = 0;i < pl.E.size();++ i){ 317 if(!P_O_IN(pl.E[i].p1,o)) return false; 318 } 319 return true; 320 }//ploygon in circle 321 322 bool O_O_IN(O o1,O o2){ 323 if(o1.r > o2.r) return false; 324 else return o2.r-o1.r <= P_P_DIST(o1.o,o2.o); 325 }//circle in circle 326 327 //================================= node ================================== 328 void P_O_ND(V p,O o,double &x,double &y){ 329 if(p == o.o){ 330 x = y = oo; 331 return; 332 } 333 L a; 334 double lth = P_O_DIST(p,o); 335 a.p1 = o.o,a.p2 = p; 336 a.getk(); 337 if(a.pa == 1){ 338 if(p.py > o.o.py) x = p.px,y = o.o.py + o.r; 339 else x = p.px,y = o.o.py - o.r; 340 } 341 else if(a.pa == 2){ 342 if(p.px > o.o.px) y = p.py,x = o.o.px + o.r; 343 else y = p.py,x = o.o.px - o.r; 344 } 345 else{ 346 double xp = o.r / sqrt(a.k*a.k + 1),yp = abs(a.k * xp); 347 x = o.o.px,y = o.o.py; 348 if(p.px > o.o.px) x += xp; 349 else x -= xp; 350 if(p.py > o.o.py) y += yp; 351 else y -= yp; 352 } 353 }//closest node of point and circle 354 355 void sameL_S_S_ND(L a,L b,double &x,double &y){ 356 V h = a.p1,v1,v2,v3,v4; 357 if(a.pa == 1) h.px += 100; 358 else if(a.pa == 2) h.py += 100; 359 else a.getk(),h.px += 100,h.py += 100 * a.k; 360 if(a.p1.MO2() > a.p2.MO2()) swap(a.p1,a.p2); 361 if(b.p1.MO2() > b.p2.MO2()) swap(b.p1,b.p2); 362 if(a.p2.MO2() > b.p2.MO2()) swap(a,b); 363 v1 = a.p1 - h,v2 = a.p2 - h,v3 = b.p1 - h,v4 = b.p2 - h; 364 double s1 = abs(CP(v1,v2)),s2 = abs(CP(v3,v4)),s3 = abs(CP(v1,v4)); 365 if(s1 + s2 < s3) x = y = - oo;//no node 366 else if(s1 + s2 > s3) x = y = oo;//oo nodes 367 else x = a.p2.px,y = a.p2.py; 368 }//node of segments which are on the same beeline 369 370 void S_SorB_ND(L a,L b,double &x,double &y){ 371 if(!((b.p0 == 1 && S_L_CR(a,b))||(b.p0 == 2 && S_S_CR(a,b)))){ 372 x = y = - oo;//no node 373 return; 374 } 375 if(a.pa == 1){ 376 if(b.pa == 1){ 377 if(b.p0 == 1) x = y = oo;//oo nodes 378 else sameL_S_S_ND(a,b,x,y); 379 } 380 else{ 381 x = a.p1.px; 382 b.getk(); 383 y = b.k * (x - b.p1.px) + b.p1.py; 384 } 385 } 386 else if(b.pa == 1){ 387 x = b.p1.px; 388 a.getk(); 389 y = a.k * (x - a.p1.px) + a.p1.py; 390 } 391 else if(a.pa == 2){ 392 if(b.pa == 2){ 393 if(b.p0 == 1) x = y = oo;//oo nodes 394 else sameL_S_S_ND(a,b,x,y); 395 } 396 else{ 397 y = a.p1.py; 398 b.getk(); 399 x = 1.0/a.k * (y - b.p1.py) + b.p1.px; 400 } 401 } 402 else if(b.pa == 2){ 403 y = b.p1.py; 404 a.getk(); 405 x = 1.0/b.k * (y - a.p1.py) + a.p1.px; 406 } 407 else{ 408 a.getk(),b.getk(); 409 if(a.k == b.k) sameL_S_S_ND(a,b,x,y); 410 else{ 411 x = (a.k*a.p1.px - b.k*b.p1.px + b.p1.py - a.p1.px) / (a.k - b.k); 412 y = a.k * (x - a.p1.px) + a.p1.py; 413 } 414 } 415 }//node of segment and segment(or beeline) 416 417 void SorB_BL_ND(L l,BL bl,V v[]){ 418 int j = 1; 419 for(int i = 0;i < bl.E.size();++ i){ 420 S_SorB_ND(l,bl.E[i],v[j].px,v[j].py); 421 if(v[j].px != oo && v[j].px != -oo) 422 j ++; 423 } 424 }//node(s) of segment(or beeline) and broken line 425 426 void SorB_PL_ND(L l,PL pl,V v[]){ 427 int j = 1; 428 for(int i = 0;i < pl.E.size();++ i){ 429 S_SorB_ND(l,pl.E[i],v[j].px,v[j].py); 430 if(v[j].px != oo && v[j].px != -oo) 431 j ++; 432 } 433 }//node(s) of segment(or beeline) and ploygon 434 435 void SorB_O_ND(L l,O o,V v[]){ 436 if(l.pa == 2 && P_O_IN(l.p1,o) && P_O_IN(l.p2,o)) return; 437 //L l0 = l;l0.pa = 1; 438 double dis = P_L_DIST(o.o,l); 439 V v1,v2; 440 if(dis > o.r) return; 441 if(l.p0 == 1){ 442 double h = sqrt(o.r * o.r - dis * dis); 443 if(h == 0){ 444 v[1].px = l.p1.px,v[1].py = o.o.py; 445 } 446 else{ 447 v1.px = v2.px = l.p1.px; 448 v1.py = v2.py = o.o.py; 449 v1.py += h,v2.py -= h; 450 if(l.pa == 1){ 451 v[1] = v1,v[2] = v2; 452 } 453 else{ 454 if((P_O_IN(l.p1,o) && P_O_IN(l.p2,o))||(!P_O_IN(l.p1,o) && !P_O_IN(l.p2,o))) return; 455 if(P_L_ON(v1,l) && P_L_ON(v2,l)){ 456 v[1] = v1,v[2] = v2; 457 } 458 else if(P_L_ON(v1,l)){ 459 v[1] = v1; 460 } 461 else if(P_L_ON(v2,l)){ 462 v[1] = v2; 463 } 464 } 465 } 466 } 467 else if(l.p0 == 2){ 468 double w = sqrt(o.r * o.r - dis * dis); 469 if(w == 0){ 470 v[1].py = l.p1.py,v[1].px = o.o.px; 471 } 472 else{ 473 v1.py = v2.py = l.p1.py; 474 v1.px = v2.px = o.o.px; 475 v1.px += w,v2.px -= w; 476 if(l.pa == 1){ 477 v[1] = v1,v[2] = v2; 478 } 479 else{ 480 if((P_O_IN(l.p1,o) && P_O_IN(l.p2,o))||(!P_O_IN(l.p1,o) && !P_O_IN(l.p2,o))) return; 481 if(P_L_ON(v1,l) && P_L_ON(v2,l)){ 482 v[1] = v1,v[2] = v2; 483 } 484 else if(P_L_ON(v1,l)){ 485 v[1] = v1; 486 } 487 else if(P_L_ON(v2,l)){ 488 v[1] = v2; 489 } 490 } 491 } 492 } 493 else{ 494 l.getk(); 495 double A = l.k * l.k + 1,B = 2 * (-l.k * l.p1.px + l.p1.py - o.o.py - o.o.px),C = (-l.k * l.p1.px + l.p1.py - o.o.py) * (-l.k * l.p1.px + l.p1.py - o.o.py) - o.r * o.r,D;// = 4 * ((-l.k * l.p1.px + l.p1.py - o.o.py - o.o.px) * (-l.k * l.p1.px + l.p1.py - o.o.py - o.o.px) - (l.k * l.k + 1) * ((-l.k * l.p1.px + l.p1.py - o.o.py) * (-l.k * l.p1.px + l.p1.py - o.o.py) - o.r * o.r)); 496 D = B * B - 4 * A * C; 497 v1.px = (-B + sqrt(D))/(2*A),v1.py = l.k * (v1.px - l.p1.px) + l.p1.py; 498 v2.px = (-B - sqrt(D))/(2*A),v2.py = l.k * (v2.px - l.p1.px) + l.p1.py; 499 if(D == 0){ 500 v[1] = v1; 501 } 502 else{ 503 if(l.pa == 1){ 504 v[1] = v1,v[2] = v2; 505 } 506 else{ 507 if((P_O_IN(l.p1,o) && P_O_IN(l.p2,o))||(!P_O_IN(l.p1,o) && !P_O_IN(l.p2,o))) return; 508 if(P_L_ON(v1,l) && P_L_ON(v2,l)){ 509 v[1] = v1,v[2] = v2; 510 } 511 else if(P_L_ON(v1,l)){ 512 v[1] = v1; 513 } 514 else if(P_L_ON(v2,l)){ 515 v[1] = v2; 516 } 517 } 518 } 519 } 520 }//node(s) of segment(or beeline) and circle 521 522 //============================ convex closure ============================= 523 bool cmp(V a,V b){ 524 return CP(a,b) > 0; 525 } 526 527 stack<V,vector<V> > GCC(V src[],int n){ 528 V p0 = src[1]; 529 int o0 = 1; 530 for(int i = 2;i <= n;++ i){ 531 if(p0.py != src[i].py){ 532 p0 = p0.py < src[i].py? p0:(o0 = i,src[i]); 533 } 534 else p0 = p0.px < src[i].py? p0:(o0 = i,src[i]); 535 } 536 V tmp[maxn]; 537 int j = 0; 538 for(int i = 1;i <= n;++ i){ 539 if(i != o0){ 540 tmp[++j] = src[i]; 541 } 542 } 543 sort(tmp+1,tmp+n,cmp); 544 stack<V,vector<V> > S; 545 S.push(p0),S.push(tmp[1]),S.push(tmp[2]); 546 V ut = tmp[1];//ut: undertop 547 for(int i = 3;i < n;++ i){ 548 while(CP(S.top() - ut,tmp[i] - S.top()) <= 0){ 549 S.pop(),S.pop(); 550 V tp = ut;ut = S.top(); 551 S.push(tp); 552 } 553 S.push(tmp[i]); 554 } 555 return S; 556 }//get convex closure 557 //================================= main ================================== 558 const double OO = 1e-8; 559 const double oo = 1e9; 560 const int maxn = 10000; 561 562 563 564 void init(){ 565 } 566 567 void work(){ 568 } 569 570 int main(){ 571 init(); 572 work(); 573 return 0; 574 }
