已解决问题:
判断点是否在线段上
判断两线段是否相交
判断点是否在多边形内
判断线段、折线、多边形是否在多变形内
判断上述是否在圆内
计算上述与线段及直线的交点
凸包
待解决:
半平面交
1 #include2 #include 3 #include 4 #include 5 #include 6 #include 7 #include 8 #include 9 #include 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 line102 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_S115 return exp <= 0;116 break;117 case 2://S_S118 return exp < 0;119 break;120 case 3://R_S121 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 cross135 };//line136 137 struct BL{138 vector E;139 };//broken line140 141 struct PL{142 vector E;143 };//ploygon144 145 struct O{146 V o;147 double r;148 };//circle149 150 //================================ data ===================================151 double CP(V a,V b){152 return a.px*b.py - a.py*b.px;153 }//cross product154 155 double DP(V a,V b){156 return a.px*b.px + a.py*b.py;157 }//dot product158 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 points162 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 line181 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 ploygon189 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 circle193 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 area203 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' cross231 232 bool S_L_CR(L a,L b){233 return b.CR(a);234 }//segment and beeline's cross235 236 bool S_R_CR(L a,L b){237 return b.CR(a);238 }//segment and ray's cross239 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 ploygon257 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;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 ploygon277 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 ploygon284 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 ploygon291 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 ploygon299 300 bool P_O_IN(V p,O o){301 return (p-o.o).MO2() <= o.r * o.r;302 }//point in circle303 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 circle307 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 circle314 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 circle321 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 circle326 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 circle354 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 node366 else if(s1 + s2 > s3) x = y = oo;//oo nodes367 else x = a.p2.px,y = a.p2.py;368 }//node of segments which are on the same beeline369 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 node373 return;374 }375 if(a.pa == 1){376 if(b.pa == 1){377 if(b.p0 == 1) x = y = oo;//oo nodes378 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 nodes394 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 line425 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 ploygon434 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 circle521 522 //============================ convex closure =============================523 bool cmp(V a,V b){524 return CP(a,b) > 0;525 }526 527 stack > 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 > S;545 S.push(p0),S.push(tmp[1]),S.push(tmp[2]);546 V ut = tmp[1];//ut: undertop547 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 closure557 //================================= 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 }