计算几何模板

来源:互联网 发布:sift算法matlab代码 编辑:程序博客网 时间:2024/05/18 12:05

废话不多说,先附上一个粗糙的板吧 ,然后慢慢完善~

const double pi=acos(-1.0);const double eps=1e-8;/*my_template*//*Computational Geometry*/const int N = 10010;struct Point{double x,y;Point(double _x=0,double _y=0){x=_x,y=_y;}void read(){scanf("%lf%lf",&x,&y);}}p[N];struct Line{double a,b,c;//表示一条直线ax+by+c=0;Line(double _a,double _b,double _c):a(_a),b(_b),c(_c){}Line(){}void read(){scanf("%lf%lf%lf",&a,&b,&c);}};struct Lineseg{Point s,e;Lineseg(){}Lineseg(Point _s ,Point _e ){s = _s; e = _e;}void read(){s.read();e.read();}};class CG2d{public:bool equal_double(double a,double b){return (fabs(a-b)<eps);}double dist_2(Point p1,Point p2){return ((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));}double dist(Point p1,Point p2){return sqrt(dist_2(p1,p2));}bool equal_point(Point p1,Point p2){return equal_double(p1.x,p2.x)&&equal_double(p1.y,p2.y);}/*Cross product ;>0 ep在矢量opsp的逆时针方向 ;=0 共线*/double multiply(Point op,Point sp ,Point ep){return ((sp.x-op.x)*(ep.y-op.y)-(ep.x-op.x)*(sp.y-op.y));}/*Dot product ;r<0 夹角为锐角; =0 直角 ;>0 钝角*/double dotmultiply(Point op ,Point sp, Point ep){return ((sp.x - op.x)*(ep.x-op.x)+(sp.y-op.y)*(ep.y-op.y));}bool onlineseg(Lineseg l,Point p){return ((multiply(l.s,l.e,p)==0)&&( ((p.x-l.s.x)*(p.x-l.e.x)<=0) && ((p.y-l.s.y)*(p.y-l.e.y)<=0) ) );}/*返回点p以点o为圆心逆时针旋转alpha(单位:弧度)后所在的位置*/Point rotate(Point o ,double alpha,Point p){Point ret;p.x -= o.x;p.y -= o.y;ret.x = p.x*cos(alpha) - p.y*sin(alpha) + o.x;ret.y = p.y*cos(alpha) + p.x*sin(alpha) + o.y;return ret;}/*angle计算夹角*/double angle(Point op,Point sp,Point ep){double cosfi,phi,norm;double dsx = sp.x - op.x,dsy = sp.y - op.y;double dex = ep.x - op.x,dey = ep.y - op.y;cosfi = dsx*dex + dsy * dey;norm = (dsx*dsx+dsy*dsy)*(dex*dex+dey*dey);cosfi /= sqrt(norm);if (cosfi >= 1.0) return 0;if (cosfi <= -1.0) return -pi;phi = acos(cosfi);if (dsx*dey-dsy*dex>0) return phi;return -phi;}/*******************线段与直线的基本运算***************//*判断点与线段的关系,用途很广泛*/  double relation(Point p,Lineseg l){return dotmultiply(l.s,p,l.e)/dist_2(l.s,l.e);}/*垂足点p到线段l的垂足*/Point perpendicular(Point p,Lineseg l){double r = relation(p,l);Point tp(l.s.x + r * (l.e.x - l.s.x),l.s.y + r * (l.e.y - l.s.y));return tp;}/*求点P到线段l的最短距离,并返回线段上距该点最近的点np*/double ptolinesegdist(Point p,Lineseg l,Point &np){double r = relation(p,l);if (r<0) np = l.s;else if (r>1) np = l.e ;else np = perpendicular(p,l);return dist(p,np);}/*求点p到线段l所在直线的距离*/double ptolinedist(Point p,Lineseg l){return fabs(multiply(l.s,l.e,p))/dist(l.s,l.e);}/*求点到折线集的最近距离,并返回最近点*/double ptopointset(int vcnt,Point pointset[],Point p,Point &q){double rd = 1000000000.0,td;   //rd初始化 赋值成一个很大的值Lineseg l;Point rq,tq;for(int i = 0;i < vcnt-1; ++i){l.s = pointset[i];l.e = pointset[i+1];td = ptolinesegdist(p,l,tq);if (td < rd) rd = td,rq = tq;}q = rq;return rd; }double ptopolygon(int vcnt,Point pointset[],Point p,Point &q){Point tq;double d = ptopointset(vcnt,pointset,p,q);double td = ptolinesegdist(p,Lineseg(pointset[vcnt-1],pointset[0]),tq);if (td < d) {q = tq;return td;}return d;}/*判断圆是否在多边形内 这里木有判断圆心是否在多边形内,默认在里面*/bool CircleInsidePolygon(int vcnt,Point polygon[],Point center,double radius){Point q;double d = ptopolygon(vcnt,polygon,center,q);if (radius <= d || fabs(d-radius) < eps) return 1;return 0;}/*返回 两个矢量l1和l2的夹角的余弦 (-1..1) 注意如果要求余弦的话,反余弦函数的定义是从0到pi的*/double cosine(Lineseg l1,Lineseg l2){return ((l1.e.x-l1.s.x) * (l2.e.x - l2.s.x) + (l1.e.y - l1.s.y) * (l2.e.y - l2.s.y))/(dist(l1.s,l1.e)*dist(l2.s,l2.e));}/*返回线段l1 与 l2之间的夹角, 单位为弧度 范围(-pi,pi)*/double lsangle(Lineseg l1,Lineseg l2){Point o(0,0),s,e;s.x = l1.e.x - l1.s.x;s.y = l1.e.y - l1.s.y;e.x = l2.e.x - l2.s.x;e.y = l2.e.y - l2.s.y;return angle(o,s,e);}/*如果线段u和v 相交(包括相交在端点处)时,然后true*/bool Isintersect(Point a1, Point a2, Point b1, Point b2) {//判断两条线段是否相交     return  min(a1.x, a2.x) <= max(b1.x, b2.x)&&         min(a1.y, a2.y) <= max(b1.y, b2.y) &&         min(b1.x, b2.x) <= max(a1.x, a2.x) &&         min(b1.y, b2.y) <= max(a1.y, a2.y) &&         multiply(a1, a2, b1) * multiply(a1, a2, b2) <= 0 &&         multiply(b1, b2, a1) * multiply(b1, b2, a2) <= 0 ;        }/*如果线段u和v 相交(不包括相交在端点处)时,然后true*/bool IsintersectNotOnline(Point a1, Point a2, Point b1, Point b2){return Isintersect(a1,a2,b1,b2) && (!onlineseg(Lineseg(b1,b2),a1)) &&(!onlineseg(Lineseg(b1,b2),a2)) && (!onlineseg(Lineseg(a1,a2),b1)) &&(!onlineseg(Lineseg(a1,a2),b2));}/*直线与线段相交模板:点(P1,P2)是否与过点(Q1,Q2)的直线相交 没有判规范 */bool intersect_line_seg(Point P1,Point P2,Point Q1,Point Q2) {   return multiply(Q1,P1,Q2)*multiply(Q1,Q2,P2)>=0;} //已知直线过的两点 求过这两点的直线方程 ax+by+c=0; a>=0;Line makeline(Point p1,Point p2){Line tl;tl.a = p2.y - p1.y;tl.b = p1.x - p2.x;tl.c = p1.y * p2.x - p1.x * p2.y;//if (tl.a<0) tl.a = -tl.a,tl.b=-tl.b,tl.c=-tl.c;return tl;}//计算直线的倾斜角double slope(Line l){if (fabs(l.a) < eps) return 0.0;if (fabs(l.b) < eps) return 1e20;  //斜率无穷大return -(l.a/l.b);}//返回直线的倾斜角double alpha(Line l){if (fabs(l.a)<eps) return 0.0;if (fabs(l.b)<eps) return pi/2;double k = slope(l);if (k>0) return atan(k);return pi+atan(k);}//求点p关于直线l 的对称点Point symmetry(Line l ,Point p){Point rp;rp.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);rp.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);return rp;}/*如果两条直线l1 和 l2 相交返回true 和交点p*/bool lineintersect(Line l1,Line l2,Point &p){double d = l1.a * l2.b - l2.a*l1.b;if (fabs(d)<eps) return 0;p.y = (l1.c * l2.a - l1.a * l2.c)/d;p.x = (l2.c * l1.b - l1.c * l2.b)/d;return 1;}//判断直线重合bool Line_Coincide(const Line &l1,const Line &l2){return equal_double(l1.a*l2.c,l2.a*l1.c) && equal_double(l1.b*l2.c,l2.b*l1.c);}/*判断点是否在多边形内 如果点在多边形上也视为在里面*/bool point_in_polygon(int vcnt,Point pointset[],Point p){int cnt = 0;Lineseg l(p,Point(1e10,p.y));Point pn = pointset[vcnt];pointset[vcnt] = pointset[0];fr(i,0,vcnt){if (onlineseg(Lineseg(pointset[i],pointset[i+1]),p)) return 1;if (equal_double(pointset[i].y ,pointset[i+1].y)) continue;int tmp = -1;if (onlineseg(l,pointset[i])) tmp=i;else if (onlineseg(l,pointset[i+1])) tmp =i+1;if (tmp != -1 && equal_double(pointset[tmp].y,min(pointset[i].y,pointset[i+1].y))) cnt++;if (tmp == -1 && Isintersect(l.s,l.e,pointset[i],pointset[i+1])) cnt++;}pointset[vcnt] = pn;if (cnt&1) return 1;return 0;}//已知三角形三个顶点的坐标,求三角形的面积     double triangleArea(Point p0, Point p1, Point p2) {  double k = p0.x * p1.y + p1.x * p2.y   + p2.x * p0.y - p1.x * p0.y - p2.x * p1.y - p0.x * p2.y;     return fabs(k/2);}  //凸包的极角序/*static bool gramhap_cmp(Point a,Point b){double k = cg2d.multiply(p[0],a,b);if (fabs(k)<eps) return cg2d.dist(a,p[0])<cg2d.dist(p[0],b);else if (k>0) return 1;else return 0;}//n为原来多边形的点的个数 返回凸包上的点的个数int graham(int n,Point chs[]){int mny = p[0].y ,idx = 0;fr(i,1,n){if (p[i].y < mny) mny = p[i].y,idx = i;else if (p[i].y == mny &&p[i].x <p[idx].x) idx = i;}Point tmp = p[idx];p[idx] = p[0];p[0] = tmp;sort(p+1,p+n,gramhap_cmp);chs[0] = p[n-1];chs[1] = p[0];int stop = 1,pos =1;while(pos<=n-1){double d = multiply(chs[stop],chs[stop-1],p[pos]);if (d<=0) chs[++stop] = p[pos++];else stop--;}return stop+1;}//计算凸包的周长double convex_circumference(int n,Point p[]){double ret = 0;fr(i,0,n-1) ret += dist(p[i],p[i+1]);ret += dist(p[n-1],p[0]);return ret;}double convex_area(int n,Point chs[]){double ret = 0;Point s = chs[0];fr(i , 1,n-1)  ret += multiply(s,chs[i],chs[i+1]);return fabs(ret)/2;}*/}cg2d;




利用这个板 可以轻松的切掉一般的计算几何题目

poj 1066 Treasure Hunt

枚举所有先线段的短点,与目标点判一次相交即可

#define sfint(x) scanf("%d",&x)#define sfint2(x,y) scanf("%d%d",&x,&y)#define sfint3(x,y,z) scanf("%d%d%d",&x,&y,&z)#define sfstr(c) scanf("%s",c)#define pfint(x) printf("%d\n",x)#define fr(i,s,n) for(int i=s;i<n;++i)#define _fr(i,n,s) for(int i=n-1;i>=s;--i)#define cl(a) memset(a,0,sizeof(a))Lineseg l[31];Point st;int n;int main(){fi;sfint(n);fr(i,0,n) l[i].read();st.read();int ans = 100;fr(i,0,n){int tmp = 0;fr(j,0,n) if (cg2d.Isintersect(st,l[i].s,l[j].s,l[j].e)) tmp++;ans =min(tmp,ans);tmp = 0;fr(j,0,n) if (cg2d.Isintersect(st,l[i].e,l[j].s,l[j].e)) tmp++;ans =min(tmp,ans);}if (n == 0) ans = 1;printf("Number of doors = %d\n",ans);}



poj 2653 Pick-up sticks

struct S{Lineseg l;int pos;}l[100010];list<S> lis;int n;void solve(){lis.clear();list<S>::iterator lit,tmpit;;lis.insert(lis.begin(),l[1]);fr(i ,2,n+1){for(lit = lis.begin();lit!=lis.end();){tmpit = lit;++lit;if(cg2d.Isintersect(tmpit->l.s,tmpit->l.e,l[i].l.s,l[i].l.e)){lis.erase(tmpit);}}lis.push_back(l[i]);}printf("Top sticks: ");for (lit = lis.begin();lit!=lis.end();){printf("%d",lit->pos);++lit;if (lit == lis.end()){printf(".\n");}else{printf(", ");}}}void inp(){fr(i,1,n+1) l[i].l.read(),l[i].pos = i;}int main(){#ifdef localfi;#endifwhile(sfint(n),n){inp();solve();}

poj1410 Intersection

Point rec[4];int main(){#ifdef localfi;#endifint t;Lineseg l;Point p,q;sfint(t);while(t--){l.read();rec[0].read();rec[2].read();double lx = min(rec[0].x,rec[2].x),hx = max(rec[0].x,rec[2].x);double ly = min(rec[0].y,rec[2].y),hy = max(rec[0].y,rec[2].y);rec[1].x = rec[0].x;rec[1].y = rec[2].y;rec[3].x = rec[2].x;rec[3].y = rec[0].y;bool flag = 0;if((l.s.x <= hx && l.s.x>= lx && l.s.y <= hy && l.s.y >= ly) &&(l.e.x <= hx && l.e.x>= lx && l.e.y <= hy && l.e.y >= ly) )flag = 1;if (!flag)fr(i,0,4){if (cg2d.Isintersect(rec[i],rec[(i+1)%4],l.s,l.e)){flag = 1;break;}}if (flag) puts("T");else puts("F");}return 0;}


POJ 1584 A Round Peg in a Ground Hole

Point center,p[100000],chs[100000];;int n;double r;bool check(){int stop =1,pos = 1;chs[0] = p[n-1],chs[1] = p[0];chs[n] = p[0];p[n] = p[0];bool flag = 0;while(pos<=n){double d = cg2d.multiply(chs[stop],chs[stop-1],p[pos]);if ((d<=0)) {chs[++stop] = p[pos++];}else {flag = 1;break;}}if (flag == 1){pos = 1;stop =1;while(pos<=n){double d = cg2d.multiply(chs[stop],chs[stop-1],p[pos]);if (d>=0) {chs[++stop] = p[pos++];}else return 0;}}return 1;}int main(){fi;while(1){sfint(n);if (n<3) break;sfdl(r);center.read();fr(i,0,n) p[i].read();if (check()){if (cg2d.point_in_polygon(n,p,center)){if (cg2d.CircleInsidePolygon(n,p,center,r)) puts("PEG WILL FIT");else puts("PEG WILL NOT FIT");}else puts("PEG WILL NOT FIT");}else{puts("HOLE IS ILL-FORMED");}}}


POJ 1113 Wall


Point chs[N],p[N];int main(){int n;double l;sfint(n);sfdl(l);fr(i , 0,n){p[i].read();}int vcnt = cg2d.graham(n,chs);double ans = cg2d.convex_circumference(vcnt,chs);ans += pi*2*l;printf("%.0f\n",ans);}


POJ 1228 Grandpa's Estate

Point chs[N];int n;void solve(){int vcnt = cg2d.graham(n,chs);bool flag = 1;Point prepre=chs[1],pre=chs[2];chs[vcnt] = chs[0];int ct = 2;fr(i , 3,vcnt+1){  //最后一个还要和起点检查一下下;if (ct == 2){if (!(fabs(cg2d.multiply(prepre,pre,chs[i]))<eps)) {flag = 0 ; break;}else {ct++;prepre = pre;pre=chs[i];}}else{if (fabs(cg2d.multiply(prepre,pre,chs[i]))<eps) ct++;else ct = 2;prepre = pre;pre = chs[i];}}if (!flag) puts("NO");else{if (fabs(cg2d.multiply(prepre,pre,chs[1]))<eps) puts("YES");else puts("NO");}}int main(){int t;sfint(t);while(t--){sfint(n);fr(i , 0,n) p[i].read();if (n<6) puts("NO");else solve();}return 0;}


POJ 3608 Bridge Across Islands

int n,m,cntp ,cntq;Point  chsp[N];Point chsq[N];double solve(Point p[],Point q[]){int sp=0,sq=0;fr(i , 1,cntp)  if ((p[i].y<p[sp].y)||(cg2d.equal_double(p[i].y,p[sp].y)&&p[i].x<p[sp].x))sp = i;fr(i , 1,cntq)if ((q[i].y>q[sq].y)||(cg2d.equal_double(q[i].y,q[sq].y)&&q[i].x>q[sq].x))sq = i;int tp = sp,tq = sq;double ans = cg2d.dist(p[sp],q[sq]);Point xxx;do{double k = cg2d.multiply(Point(0,0),Point(p[(sp+1)%cntp].x-p[sp].x,p[(sp+1)%cntp].y-p[sp].y),Point(q[(sq+1)%cntq].x-q[sq].x,q[(sq+1)%cntq].y-q[sq].y));if (cg2d.equal_double(k,0)){ans = min(ans,cg2d.ptolinesegdist(p[sp],Lineseg(q[sq],q[(sq+1)%cntq]),xxx));ans = min(ans,cg2d.ptolinesegdist(p[(sp+1)%cntp],Lineseg(q[sq],q[(sq+1)%cntq]),xxx));ans = min(ans,cg2d.ptolinesegdist(q[sq],Lineseg(p[sp],p[(sp+1)%cntp]),xxx));ans = min(ans,cg2d.ptolinesegdist(q[(sq+1)%cntq],Lineseg(p[sp],p[(sp+1)%cntp]),xxx));sp = (sp+1)%cntp;sq = (sq+1)%cntq;}else if (k<0){ans = min(ans,cg2d.ptolinesegdist(q[sq],Lineseg(p[sp],p[(sp+1)%cntp]),xxx));sp = (sp+1)%cntp;}else{ans = min(ans,cg2d.ptolinesegdist(p[sp],Lineseg(q[sq],q[(sq+1)%cntq]),xxx));sq = (sq+1)%cntq;}}while(!(tp == sp&&tq == sq));return ans;}int main(){while(1){sfint2(n,m);if (n==0 && m==0) break;fr(i, 0,n) p[i].read();fr(i,0,m) q[i].read();cntp = cg2d.graham(n,p,chsp);//求最远点对时graham:pos<n-1cntq = cg2d.graham(m,q,chsq);printf("%.5f\n",solve(chsp,chsq));}return 0;}

原创粉丝点击