Computational Geometry Templates
2016-02-14 22:08
429 查看
//author: birdstorm #include<bits/stdc++.h> using namespace std; /*************useful TIPS in complex*********//** *1. use eps carefully *2. remember useful functions in complex: * abs() for distance, * arg() for angle, * conj() for conjunction, * norm() for |v|^2, * polar(double mag, double ang=0.0) for creating complex * real(), imag(), ... *3. useful functions 2: * cos(), sin(), acos(), asin(), tan(), atan(), cosh(), sinh(), tanh(), * pow(), exp(), sqrt(), log(), ... * can also be used in complex directly! *4. debug with cout<<v<<endl; *5. we can also use "type _Complex" in definition(C99/C++11) * e.g. double _Complex p; float _Complex q; *6. To Be Continued... **************************************************/ /*********CONSTANTS**************/ const double eps=1.0e-8; const double pi=acos(-1.0); /***************DEFINATION*********************/ typedef complex<double> point; typedef vector<point> polygon; typedef point Vec; struct line : public vector<point> { line() {} line(const point& a, const point& b) { push_back(a); push_back(b); } }; typedef vector<line> Obstacle; struct Circle{ point center; double radius; Circle(const point& c, const double& r):center(c),radius(r) {} }; struct polygon_convex{ polygon p; polygon_convex(int sz=0){ p.resize(sz); } }; /******************COMMONLY-USED FUNCTIONS*************/ point unit(const point& v) { return v/abs(v);} point ortho(const point& v) { return v*point(0,1);} point spin(const point& v, const double &angle) {return v*point(cos(angle),sin(angle));} inline point vec(const line& l) { return l[1]-l[0];} bool equal(const double a, const double b) { return abs(a-b)<eps;} bool equal(const point& a, const point& b) { return abs(a-b)<eps;} inline double dot (const point& a, const point& b) { return (a*conj(b)).real();} inline double cross(const point& a, const point& b) { return (conj(a)*b).imag();} inline int dlcmp(const double& x){ return x<-eps?-1:x>eps; } /**************COMMONLY-USED FUNCTIONS 2***************/ inline int ccw(const point& a, const point& b, const point& c) { point u=b-a, v=c-a; if(cross(u,v)>0) return 1; if(cross(u,v)<0) return -1; if(dot(u,v)<0) return 2; if(abs(u)<abs(v)) return -2; return 0; } inline int ccw(const line& s, const point& p) { return ccw(s[0], s[1], p); } /**************COMPARISON********************/ namespace std{ bool operator<(const point& a, const point& b) { return a.real()!=b.real()?a.real()<b.real():a.imag()<b.imag(); } bool operator==(const point& a, const point& b) { return abs(a-b)<eps; } } bool cmp_less(const point& a, const point& b){ return a.real()!=b.real()?a.real()<b.real():a.imag()<b.imag(); } /**************INTERSECTIONS Judging****************/ bool intersectLP(const line& l, const point& p) { return abs(cross(l[1]-p, l[0]-p)) < eps; } bool intersectSP(const line& s, const point& p) { return abs(s[0]-p)+abs(s[1]-p) < abs(s[1]-s[0])+eps; } bool intersectLL(const line& l, const line& m) { return cross(vec(l), vec(m)) > eps || cross(vec(l), m[0]-l[0]) < eps; } bool intersectLS(const line& l, const line& s) { return cross(vec(l),s[0]-l[0]) * cross(vec(l),s[1]-l[0]) <= 0; } bool intersectSS(const line& s, const line& t) { return ccw(s,t[0])*ccw(s,t[1]) <= 0 && ccw(t,s[0])*ccw(t,s[1]) <= 0; } /************INTERSECTIONS Calculating***************/ point crosspoint(const line& l, const line& m) { ///use Intersection Judging first double A = cross(vec(l), vec(m)); double B = cross(vec(l), l[1]-m[0]); if(abs(A)<eps) { // parallel return m[0]; // sameline } return m[0] + B/A*vec(m); } point intersection_point(const line& A, const Circle& C){ double dx=real(A[1]-A[0]), dy=imag(A[1]-A[0]); double rx=real(A[0]-C.center), ry=imag(A[0]-C.center); double a=dx*dx+dy*dy, b=2*dx*rx+2*dy*ry, c=rx*rx+ry*ry-C.radius*C.radius; double delta=b*b-4*a*c; if(delta<-eps) return point(-10000,-10000); double ans=(-b+sqrt(delta))/(2.0*a); return point(A[0].real()+ans*dx,A[0].imag()+ans*dy); } point projection(const line &l, const point &p) { double t = dot(p-l[0], l[0]-l[1]) / norm(l[0]-l[1]); return l[0] + t*(l[0]-l[1]); } double distanceSP(const line &s, const point &p) { const point r = projection(s, p); if (intersectSP(s, r)) return abs(r - p); return min(abs(s[0] - p), abs(s[1] - p)); } line generateSP(const line &s, const point &p) { const point r = projection(s, p); if(intersectSP(s, r)) return line(r, p); else{ if(abs(s[0]-p) < abs(s[1] - p)) return line(s[0], p); else return line(s[1], p); } return s; } bool isSameLine(const point& a, const point& b, const line& l){ double x = distanceSP(l, a); double y = distanceSP(l, b); return x < eps && y < eps; } /**************more complicated APPLICATION...***********/ bool comp(const point& a, const point& b){ return abs(a-b) < eps; } vector<point> crosspointSO(const line& sight, const vector<line>& obstacle){ vector<point> res; for(int i=0; i<obstacle.size(); i++) { if(intersectSS(sight, obstacle[i])) res.push_back(crosspoint(sight, obstacle[i])); } sort(res.begin(), res.end()); res.erase(unique(res.begin(), res.end(), comp), res.end()); if(res.size() <= 1) return vector<point>(); for(int i=0; i<obstacle.size(); i++) { if(isSameLine(res[0], res[1], obstacle[i])) return vector<point>(); } return res; } line calcSight(const line& sight, const vector<Obstacle>& obstacles){ vector<pair<double,point> > res; for(int i=0; i<obstacles.size(); i++) { vector<point> p = crosspointSO(sight, obstacles[i]); for(int j=0; j<p.size(); j++) { res.push_back(make_pair(abs(p[j]-sight[0]), p[j])); } } sort(res.begin(), res.end()); if(res.size() == 0) return sight; return line(sight[0], res[0].second); } bool reachable(const line& sight, const vector<Obstacle>& obstacles){ for(int i=0; i<obstacles.size(); i++) { vector<point> p = crosspointSO(sight, obstacles[i]); if(p.size() != 0) return false; } return true; } double TriAngleCircleInsection(Circle C, point A, point B){ Vec OX = A-C.center, OY = B-C.center; Vec YX = A-B, YO = C.center-B; Vec XY = B-A, XO = C.center-A; double da = abs(OX), db = abs(OY),dc = abs(XY), r = C.radius; if(dlcmp(cross(OX,OY)) == 0) return 0; if(dlcmp(da-r) < 0 && dlcmp(db-r) < 0) return cross(OX,OY)*0.5; else if(db < r && da >= r) { double x = (dot(YX,YO) + sqrt(r*r*dc*dc-cross(YX,YO)*cross(YX,YO)))/dc; double TS = cross(OX,OY)*0.5; return asin(TS*(1-x/dc)*2/r/da)*r*r*0.5+TS*x/dc; } else if(db >= r && da < r) { double y = (dot(XY,XO)+sqrt(r*r*dc*dc-cross(XY,XO)*cross(XY,XO)))/dc; double TS = cross(OX,OY)*0.5; return asin(TS*(1-y/dc)*2/r/db)*r*r*0.5+TS*y/dc; } else if(fabs(cross(OX,OY)) >= r*dc || dot(XY,XO) <= 0 || dot(YX,YO) <= 0) { if(dot(OX,OY) < 0) { if(cross(OX,OY) < 0) return (-acos(-1.0)-asin(cross(OX,OY)/da/db))*r*r*0.5; else return ( acos(-1.0)-asin(cross(OX,OY)/da/db))*r*r*0.5; } else return asin(cross(OX,OY)/da/db)*r*r*0.5; } else { double x = (dot(YX,YO)+sqrt(r*r*dc*dc-cross(YX,YO)*cross(YX,YO)))/dc; double y = (dot(XY,XO)+sqrt(r*r*dc*dc-cross(XY,XO)*cross(XY,XO)))/dc; double TS = cross(OX,OY)*0.5; return (asin(TS*(1-x/dc)*2/r/da)+asin(TS*(1-y/dc)*2/r/db))*r*r*0.5 + TS*((x+y)/dc-1); } } /********************POLYGONS**********************/ polygon_convex convex_hull(polygon a){ polygon_convex res(2*a.size()+5); sort(a.begin(),a.end(),cmp_less); int m=0, n=a.size(); for(int i=0; i<n; i++){ while(m>1&&dlcmp(cross(res.p[m-1]-res.p[m-2],a[i]-res.p[m-2]))<=0) --m; res.p[m++]=a[i]; } int k=m; for(int i=n-2;i>=0;--i){ while(m>k&&dlcmp(cross(res.p[m-1]-res.p[m-2],a[i]-res.p[m-2]))<=0) --m; res.p[m++]=a[i]; } res.p.resize(m); if(a.size()>1) res.p.resize(m-1); return res; } bool contain(const polygon_convex& a, const point& b){ int n=a.p.size(); #define next(i) ((i+1)%n) int sign=0; for(int i=0; i<n; i++){ int x=dlcmp(cross(a.p[i]-b,a.p[next(i)]-b)); if(x){ if(sign){ if(sign!=x) return false; } else sign=x; } } return true; } /******************END********************/ int main(){ return 0; } 圆交多边形: #include <bits/stdc++.h> using namespace std; typedef complex<double> point; const double pi=acos(-1.0); const double eps=1.0e-8; int dcmp(double k){ return k<-eps?-1:k>eps; } const int MAXN=105; double cross(const point& a, const point& b){ return (conj(a)*b).imag(); } double dot (const point& a, const point& b){ return (a*conj(b)).real(); } point crosspt(const point& a, const point& b, const point& p, const point& q){ double a1=cross(b-a,p-a); double a2=cross(b-a,q-a); return (p*a2-q*a1)/(a2-a1); } double mysqrt(double x){ return sqrt(max(0.0,x)); } double sqr(double x){ return x*x; } double sector_area(const point& a, const point& b, const double& r){ double theta=arg(a)-arg(b); while(theta<=0) theta+=2*pi; while(theta>2*pi) theta-=2*pi; theta=min(theta,2*pi-theta); return r*r*theta/2; } void circle_cross_line(point a, point b, point o, double r, point p[], int& num){ double x0=o.real(), y0=o.imag(); double x1=a.real(), y1=a.imag(); double x2=b.real(), y2=b.imag(); double dx=x2-x1, dy=y2-y1; double A=dx*dx+dy*dy; double B=2*dx*(x1-x0)+2*dy*(y1-y0); double C=sqr(x1-x0)+sqr(y1-y0)-sqr(r); double delta=B*B-4*A*C; num=0; if(dcmp(delta)>=0){ double t1=(-B-mysqrt(delta))/(2*A); double t2=(-B+mysqrt(delta))/(2*A); if(dcmp(t1-1)<=0&&dcmp(t1)>=0){ p[num++]=point(x1+t1*dx,y1+t1*dy); } if(dcmp(t2-1)<=0&&dcmp(t2)>=0){ p[num++]=point(x1+t2*dx,y1+t2*dy); } } } double calc(const point& a, const point& b, const double& r){ point p[2]; int num=0; int ina=dcmp(abs(a)-r)<0; int inb=dcmp(abs(b)-r)<0; if(ina){ if(inb){ return abs(cross(a,b))/2.0; } else{ circle_cross_line(a,b,point(0,0),r,p,num); return sector_area(b,p[0],r)+abs(cross(a,p[0]))/2.0; } } else{ if(inb){ circle_cross_line(a,b,point(0,0),r,p,num); return sector_area(p[0],a,r)+abs(cross(p[0],b))/2.0; } else{ circle_cross_line(a,b,point(0,0),r,p,num); if(num==2){ return sector_area(a,p[0],r)+sector_area(p[1],b,r)+abs(cross(p[0],p[1]))/2.0; } else{ sector_area(a,b,r); } } } } point p[MAXN]; int n; double area(const double& r){ double ret=0; for(int i=0; i<n; i++){ int sgn=dcmp(cross(p[i],p[i+1])); if(sgn){ ret+=sgn*calc(p[i],p[i+1],r); } } return ret; } int main(){ double x0, y0, v0, theta, t, g, R; while(~scanf("%lf%lf%lf%lf%lf%lf%lf",&x0,&y0,&v0,&theta,&t,&g,&R)){ if(dcmp(x0)==0&&dcmp(y0)==0&&dcmp(R)==0&&dcmp(g)==0&&dcmp(t)==0&&dcmp(theta)==0&&dcmp(v0)==0) break; theta=theta/180.0*pi; double dx=v0*cos(theta), dy=v0*sin(theta); x0+=dx*t; y0+=dy*t-0.5*g*t*t; cin>>n; double x, y; for(int i=0; i<n; i++){ scanf("%lf%lf",&x,&y); x-=x0; y-=y0; p[i]=point(x,y); } p =p[0]; printf("%.2f\n",abs(area(R))); } return 0; }
相关文章推荐
- makefile_2
- Exchange 2010 (三) HUB NLB部署
- Android SDK Manager 更新代理配置
- 查看windows系统配置的方法
- 从头认识Spring-2.3 注解装配-@autowired(5)-限定器@Qualifier(1)
- Linux中du命令:同样可以查看使用的空间,但是与df有不同
- iOS开发之网络编程--XCode7 更新以来需要手动设置的内容
- hadoop(教程,资料参考网站)
- 把REACT-NATIVE需要的安装程序都打了个包,就不用网上下了
- 线程基础:JDK1.5+(9)——线程新特性(中)
- 有关eclipse for java ee版本遇到的坑( Context initialization failed)
- poj 3132
- Multiple representations of the same entity are being merged解决方法
- 学习Discuz! X3.2记录:修改标签“Powered by Discuz!”的一种方法
- eclipse中离线安装genymotion插件
- poj 1050(矩阵求和问题dp)
- 让jsp页面自动跳转方法
- response.sendRedirect 传递参数的问题
- 文件存储的几种方式总结
- Exchange 2010 (二) MBX DAG部署