标签:
写几何题总是提心吊胆。精度问题真心吓人。
其实思路挺简单的一道题,真是什么算法和几何double搞到一块,心里就虚虚的。
思路:求出所有圆之间的交点,然后用这些交点跑一遍最短路就可以了。
Time Limit: 10000/3000 MS (Java/Others) Memory Limit: 65536/32768 K (Java/Others)
Total Submission(s): 1244 Accepted Submission(s): 304
#include <iostream> #include <cmath> #include <vector> #include <string.h> #include <stdlib.h> #include <algorithm> using namespace std; #define MAX_N 110 /*------------------常量区-------------------*/ const double INF = 1e10; // 无穷大 const double EPS = 1e-8; // 计算精度 const double PI = acos(-1.0);// PI const int PINXING = 0; // 平行 const int XIANGJIAO = 1; // 相交 const int XIANGLI = 0; // 相离 const int GONGXIAN = 2; // 共线 const int CHONGDIE = -1; // 重叠 const int INSIDE = 1; // 点在图形内部 const int OUTSIDE = 0; // 点在图形外部 const int BORDER = 2; // 点在图形边界 /*-----------------类型定义区----------------*/ struct Point { // 二维点或矢量 double x, y; //double angle, dis; Point() {} Point(double x0, double y0): x(x0), y(y0) {} void read() { scanf("%lf%lf",&x,&y); } }; struct Point3D { //三维点或矢量 double x, y, z; Point3D() {} Point3D(double x0, double y0, double z0): x(x0), y(y0), z(z0) {} }; struct Line { // 二维的直线或线段 Point p1, p2; Line() {} Line(Point p10, Point p20): p1(p10), p2(p20) {} void read() { scanf("%lf%lf%lf%lf",&p1.x,&p1.y,&p2.x,&p2.y); } }; struct Line3D { // 三维的直线或线段 Point3D p1, p2; Line3D() {} Line3D(Point3D p10, Point3D p20): p1(p10), p2(p20) {} }; struct Rect { // 用长宽表示矩形的方法 w, h分别表示宽度和高度 double w, h; Rect() {} Rect(double _w,double _h) : w(_w),h(_h) {} }; struct Rect_2 { // 表示矩形,左下角坐标是(xl, yl),右上角坐标是(xh, yh) double xl, yl, xh, yh; Rect_2() {} Rect_2(double _xl,double _yl,double _xh,double _yh) : xl(_xl),yl(_yl),xh(_xh),yh(_yh) {} }; struct Circle { //圆 Point c; double r; Circle() {} Circle(Point _c,double _r) :c(_c),r(_r) {} }; typedef vector<Point> Polygon; // 二维多边形 typedef vector<Point> Points; // 二维点集 /*-------------------基本函数区---------------------*/ inline double max(double x,double y) { return x > y ? x : y; } inline double min(double x, double y) { return x > y ? y : x; } inline bool ZERO(double x) // x == 0 { return (fabs(x) < EPS); } inline bool ZERO(Point p) // p == 0 { return (ZERO(p.x) && ZERO(p.y)); } inline bool ZERO(Point3D p) // p == 0 { return (ZERO(p.x) && ZERO(p.y) && ZERO(p.z)); } inline bool EQ(double x, double y) // eqaul, x == y { return (fabs(x - y) < EPS); } inline bool NEQ(double x, double y) // not equal, x != y { return (fabs(x - y) >= EPS); } inline bool LT(double x, double y) // less than, x < y { return ( NEQ(x, y) && (x < y) ); } inline bool GT(double x, double y) // greater than, x > y { return ( NEQ(x, y) && (x > y) ); } inline bool LEQ(double x, double y) // less equal, x <= y { return ( EQ(x, y) || (x < y) ); } inline bool GEQ(double x, double y) // greater equal, x >= y { return ( EQ(x, y) || (x > y) ); } // 输出浮点数前,防止输出-0.00调用该函数进行修正! inline double FIX(double x) { return (fabs(x) < EPS) ? 0 : x; } /*------------------二维矢量运算重载区---------------------*/ bool operator==(Point p1, Point p2) { return ( EQ(p1.x, p2.x) && EQ(p1.y, p2.y) ); } bool operator!=(Point p1, Point p2) { return ( NEQ(p1.x, p2.x) || NEQ(p1.y, p2.y) ); } bool operator<(Point p1, Point p2) { if (NEQ(p1.x, p2.x)) { return (p1.x < p2.x); } else { return (p1.y < p2.y); } } Point operator+(Point p1, Point p2) { return Point(p1.x + p2.x, p1.y + p2.y); } Point operator-(Point p1, Point p2) { return Point(p1.x - p2.x, p1.y - p2.y); } double operator*(Point p1, Point p2) // 计算叉乘 p1 × p2 { return (p1.x * p2.y - p2.x * p1.y); } double operator&(Point p1, Point p2) { // 计算点积 p1·p2 return (p1.x * p2.x + p1.y * p2.y); } double Norm(Point p) // 计算矢量p的模 { return sqrt(p.x * p.x + p.y * p.y); } /*-------------------基本函数区------------------*/ //得到两点之间的距离 double Dis(Point p1,Point p2) { return sqrt( (p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y) ); } //求二维平面上点到直线的距离 double Dis(Point p, Line L) { return ( fabs((p - L.p1) * (L.p2 - L.p1)) / Norm(L.p2 - L.p1) ); } //得到两点之间距离的平方,为减少误差用 double Dis2(Point p1,Point p2) { return (p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y); } // 把矢量p旋转角度angle (弧度表示) // angle > 0表示逆时针旋转 // angle < 0表示顺时针旋转 Point Rotate(Point p, double angle) { Point result; result.x = p.x * cos(angle) - p.y * sin(angle); result.y = p.x * sin(angle) + p.y * cos(angle); return result; } //得到向量p与x正半轴的夹角[0,2PI) double GetAngle(Point p) { double tmp=atan2(p.y,p.x); if(tmp<0) tmp=2*PI+tmp; return tmp; } //得到两个向量之间的夹角[0,PI] //若p1按顺时针转到p2的角小于PI(p1在p2的逆时针方向),则返回正数.否则返回负数. double GetAngle(Point p1,Point p2) { double tmp = GetAngle(p1) - GetAngle(p2); if( GT( tmp,PI) ) return -(2*PI-tmp); if( LEQ( tmp,-PI) ) return (tmp+2*PI); return tmp; } // 判断二维平面上点是否在线段上 // 输入:任意点p,和任意直线L // 输出:p在线段上返回1,否则返回0 bool OnSeg(Point p, Line L) { return ( ZERO( (L.p1 - p) * (L.p2 - p) ) && LEQ((p.x - L.p1.x)*(p.x - L.p2.x), 0) && LEQ((p.y - L.p1.y)*(p.y - L.p2.y), 0) ); } // 判断二维平面上点p是否在直线L上,在线段上返回1,否则返回0 bool OnLine(Point p, Line L) { return ZERO( (p - L.p1) * (L.p2 - L.p1) ); } bool OnCir(Point p,Circle cir) { return EQ( (p.x-cir.c.x)*(p.x-cir.c.x)+(p.y-cir.c.y)*(p.y-cir.c.y),cir.r*cir.r ); } //得到点p到直线L的距离,并返回p到直直线L的最近点rep double PtoLine(Point p,Line L,Point& rep) { if(L.p1==L.p2) { rep=L.p1; return Dis(p,L.p1); } Point a,b; a = L.p2-L.p1; b = p-L.p1; double dis12 = Dis(L.p1,L.p2); double dis = ( fabs(a*b) )/dis12; double k = (a&b)/(Norm(a)*dis12) ; rep.x = L.p1.x + k*(L.p2.x-L.p1.x); rep.y = L.p1.y + k*(L.p2.y-L.p1.y); return dis; } //得到点P到线段L的距离,并放回p到线段L的最近点rep double PointToSeg(Point P, Line L,Point& rep) { if(L.p1 == L.p2) { rep = L.p1; return Dis(rep,P);//如果线段是一个点,返回这个点。 } Point result; double a, b, t; a = L.p2.x - L.p1.x; b = L.p2.y - L.p1.y; t = ( (P.x - L.p1.x) * a + (P.y - L.p1.y) * b ) / (a * a + b * b);//线段上的投影 if ( GEQ(t, 0) && LEQ(t, 1) ) { result.x = L.p1.x + a * t;//值得学习的由比例求坐标的方法。 result.y = L.p1.y + b * t; } else { if ( Norm(P - L.p1) < Norm(P - L.p2) ) { result = L.p1; } else { result = L.p2; } } return Dis(result, P); } /*-------------------几何题面积计算(注意正负!)----------------------*/ // 根据三个顶点坐标计算三角形面积 // 面积的正负按照右手旋规则确定,向量AB->向量AC double Area(Point A, Point B, Point C) { return ((B-A)*(C-A) / 2.0); } // 根据三条边长计算三角形面积 double Area(double a, double b, double c) { double s = (a + b + c) / 2.0; return sqrt(s * (s - a) * (s - b) * (s - c)); } //求圆的面积 double Area(Circle C) { return PI * C.r * C.r; } // 计算多边形面积,复杂度:O(顶点数) // 面积的正负按照右手旋规则确定,顺时针为负 double Area(Polygon _poly) { int nsize=_poly.size(); double area=0; for(int i=0;i<nsize;i++) { area += _poly[i]*_poly[(i+1)%nsize]; } return area/2.0; } //求两条直线之间的关系(二维) //输入:两条不为点的直线 //输出:相交返回XIANGJIAO和交点p,平行返回PINGXING,共线返回GONGXIAN int LineAndLine(Line L1,Line L2,Point &p) { Point px,py; px = L1.p1 - L1.p2; py = L2.p1 - L2.p2; if( EQ(px*py,0) )//平行或者共线 { if( ZERO( (L2.p1-L1.p1)*py ) ) //共线 { return GONGXIAN; } return PINXING; } double xa,xb,xc,ya,yb,yc; xa=(L1.p2.y-L1.p1.y); xb=(L1.p1.x-L1.p2.x); xc=(L1.p1.y*L1.p2.x-L1.p1.x*L1.p2.y); ya=(L2.p2.y-L2.p1.y); yb=(L2.p1.x-L2.p2.x); yc=(L2.p1.y*L2.p2.x-L2.p1.x*L2.p2.y); p.y = (xa*yc-xc*ya)/(xb*ya-xa*yb); p.x = (xb*yc-xc*yb)/(xa*yb-xb*ya); return XIANGJIAO; } //判断两条线段是否相交,相交返回1 bool SegAndSeg(Line L1,Line L2) { return ( GEQ( max(L1.p1.x, L1.p2.x), min(L2.p1.x, L2.p2.x) ) && GEQ( max(L2.p1.x, L2.p2.x), min(L1.p1.x, L1.p2.x) ) && GEQ( max(L1.p1.y, L1.p2.y), min(L2.p1.y, L2.p2.y) ) && GEQ( max(L2.p1.y, L2.p2.y), min(L1.p1.y, L1.p2.y) ) && LEQ( ((L2.p1 - L1.p1) * (L1.p2 - L1.p1)) * ((L2.p2 - L1.p1) * (L1.p2 - L1.p1)), 0 ) && LEQ( ((L1.p1 - L2.p1) * (L2.p2 - L2.p1)) * ((L1.p2 - L2.p1) * (L2.p2 - L2.p1)), 0 ) ); } //求两条线段交点(二维) //输入:两条不为点的直线 //输出:相交返回XIANGJIAO和交点p,相离返回XIANGLI,重叠返回CHONGDIE int SegAndSeg(Line L1,Line L2,Point &p) { double signx,signy; //跨立实验 if( LEQ(signx=( ((L1.p2-L1.p1)*(L1.p1-L2.p1))*((L1.p2-L1.p1)*(L1.p1-L2.p2)) ),0) && LEQ(signy=( ((L2.p2-L2.p1)*(L2.p1-L1.p1))*((L2.p2-L2.p1)*(L2.p1-L1.p2)) ),0) ) { if( ZERO(signx) && ZERO(signy) ) { //线段共线 signx = min( max(L1.p1.x,L1.p2.x),max(L2.p1.x,L2.p2.x) )- max( min(L1.p1.x,L1.p2.x),min(L2.p1.x,L2.p2.x) ); signy = min( max(L1.p1.y,L1.p2.y),max(L2.p1.y,L2.p2.y) )- max( min(L1.p1.y,L1.p2.y),min(L2.p1.y,L2.p2.y) ); if( ZERO(signx) && ZERO(signy) ) //说明共线,且相交一点 { if(L1.p1==L2.p1||L1.p1==L2.p2) p=L1.p1; if(L1.p2==L2.p1||L1.p2==L2.p2) p=L1.p2; return XIANGJIAO; } else if( GEQ(signx, 0) && GEQ(signy, 0) ) { return CHONGDIE;//重叠 } else { return XIANGLI;//相离 } } return LineAndLine(L1, L2, p);//转化为直线相交 } return XIANGLI;//相离 } // 判断点p是否在简单多边形poly内, 多边形可以是凸的或凹的 // poly的顶点数目要大于等于3 // 返回值为: // INSIDE -- 点在poly内 // BORDER -- 点在poly边界上 // OUTSIDE -- 点在poly外 int InPolygon(const Polygon poly, Point p) { int i, n, count; Line ray, side; n = poly.size(); count = 0; ray.p1 = p; ray.p2.y = p.y; ray.p2.x = - INF;// 设定一个极大值 for (i = 0; i < n; i++) { side.p1 = poly[i]; side.p2 = poly[(i+1)%n]; if( OnSeg(p, side) ) { return BORDER; } // 如果side平行x轴则不作考虑 if ( EQ(side.p1.y, side.p2.y) ) { continue; } if (OnSeg(side.p1, ray)) { if ( GT(side.p1.y, side.p2.y) ) count++; } else if (OnSeg(side.p2, ray)) { if ( GT(side.p2.y, side.p1.y) ) count++; } else if ( SegAndSeg(ray, side) ) { count++; } } return ((count % 2 == 1) ? INSIDE : OUTSIDE); } //得到直线与圆的交点 int LineToCir(Line L,Circle R,Point p[2]) { if(L.p1 == L.p2)//当直线为1个点时 { if( EQ( Dis(L.p1, R.c),R.r ) ) { p[0]=L.p1; return 1; } else return 0;//相离 } Point tp;//表示圆心在直线L上的投影。 double dis=PtoLine(R.c, L,tp ); if( LT(R.r, dis) )//相离 { return 0; } if( EQ(dis,R.r) )//相切 { p[0]=tp; return 1; } double len=sqrt(R.r*R.r-dis*dis); Point onep=L.p2-L.p1; double _t = len/Norm(onep); p[0].x =tp.x + onep.x*_t; p[0].y =tp.y + onep.y*_t; onep=L.p1-L.p2; p[1].x =tp.x + onep.x*_t; p[1].y =tp.y + onep.y*_t; return 2; } //得到三角形外接圆 //注意:A,B,C三点不能共线 Circle OutCircle(Point A,Point B,Point C) { Circle tmp; double a, b, c, c1, c2; double xA, yA, xB, yB, xC, yC; a = Dis(A, B); b = Dis(B, C); c = Dis(C, A); //根据 S = a * b * c / R / 4;求半径 R tmp.r = (a*b*c)/( fabs(Area(A,B,C)) *4.0); xA = A.x; yA = A.y; xB = B.x; yB = B.y; xC = C.x; yC = C.y; c1 = (xA*xA+yA*yA - xB*xB-yB*yB) / 2; c2 = (xA*xA+yA*yA - xC*xC-yC*yC) / 2; tmp.c.x = (c1*(yA - yC)-c2*(yA - yB)) / ((xA - xB)*(yA - yC)-(xA - xC)*(yA - yB)); tmp.c.y = (c1*(xA - xC)-c2*(xA - xB)) / ((yA - yB)*(xA - xC)-(yA - yC)*(xA - xB)); return tmp; } //得到三角形内切圆 //注意:A,B,C三点不能共线 Circle InCircle(Point A,Point B,Point C) { Circle rec; double a=Dis(B,C); double b=Dis(A,C); double c=Dis(A,B); rec.c.x = (a*A.x+b*B.x+c*C.x)/(a+b+c); rec.c.y = (a*A.y+b*B.y+c*C.y)/(a+b+c); rec.r = 2*fabs( Area(A,B,C) )/(a+b+c); return rec; } //得到两圆的面积并 double CirArea(Circle _c1,Circle _c2) { if(_c2.r<_c1.r) swap(_c1,_c2); //保证_c2的半径大 double d12 = Dis(_c1.c,_c2.c); if( LEQ(_c1.r+_c2.r,d12) )//相离 { return 0; } if( LEQ(d12+_c1.r, _c2.r) )//包含 { return Area(_c1); } //相交 double _area=0; _area -= 2.0*Area(_c1.r,_c2.r,d12); double ang1 = acos( (d12*d12+_c1.r*_c1.r-_c2.r*_c2.r) / (2*d12*_c1.r) ); double ang2 = acos( (d12*d12+_c2.r*_c2.r-_c1.r*_c1.r) / (2*d12*_c2.r) ); _area += ang1*_c1.r*_c1.r+ang2*_c2.r*_c2.r; return _area; } /*---------------------代码区---------------------------*/ int CirAndCir(Circle _c1,Circle _c2,Point p[2]) { if(_c2.r < _c1.r) swap(_c1,_c2); //保证_c2的半径大 double d12 = Dis(_c1.c,_c2.c); if( LT(_c1.r+_c2.r,d12) )//相离 { return 0; } if( LT(d12+_c1.r, _c2.r) )//包含 { return 0; } if(_c1.c == _c2.c)//两个圆重叠 { return -1; } Point u,v; double t; double r1=_c1.r,r2=_c2.r; t=( 1+(r1*r1-r2*r2)/(Dis(_c1.c,_c2.c)*Dis(_c1.c,_c2.c)) ) /2; u.x = _c1.c.x + (_c2.c.x-_c1.c.x)*t; u.y = _c1.c.y + (_c2.c.y-_c1.c.y)*t; v.x = u.x + _c1.c.y - _c2.c.y; v.y = u.y - _c1.c.x + _c2.c.x; Line _l(u,v); return LineToCir(_l,_c1,p); } Circle cs[30]; Point ps[2020]; struct node { double l,r; }ls[50]; double mat[2020][2020]; int cmp(node x,node y) { if( EQ(x.l,y.l) ) return x.r<y.r; return x.l<y.l; } double Dij(double mat[2020][2020],int s,int t,int n) { double dis[2020]; int mark[2020]; for(int i=0;i<n;i++) { mark[i]=0; dis[i]=INF; } dis[s]=0; for(int ii=0;ii<n;ii++) { double mi=INF; int saveid=-1; for(int i=0;i<n;i++) { if(mark[i]==0 && dis[i] < mi ) { mi = dis[i]; saveid=i; } } if( mi > INF-1 ) { return -1; } if(saveid==t) return mi; mark[saveid] = 1; for(int i=0;i<n;i++) if( mark[i]==0 && dis[i] > dis[saveid]+mat[saveid][i] ) { dis[i] = dis[saveid]+mat[saveid][i]; } } return -1; } int main() { int T; cin>>T; int tt=1; while(T--) { int n; scanf("%d",&n); for(int i=0;i<n;i++) { scanf("%lf%lf%lf",&cs[i].c.x,&cs[i].c.y,&cs[i].r); } int id=0; ps[id++]=cs[0].c; ps[id++]=cs[n-1].c; for(int i=0;i<n;i++) { for(int j=i+1;j<n;j++) { Point p[2]; int cnt=CirAndCir(cs[i],cs[j],p); if(cnt>0)//有交点 { for(int k=0;k<cnt;k++) ps[ id++ ] = p[k]; } } } //所有圆的交点求好了 for(int i=0;i<id;i++) for(int j=0;j<id;j++) mat[i][j] = INF; for(int i=0;i<id;i++) { for(int j=i+1;j<id;j++) { if(ps[i]==ps[j]) { //两个相等的话 mat[i][j]=mat[j][i]=0; continue; } Line l(ps[i],ps[j]); int lcnt=0; //关键的一步,这条线段是否被圆全部覆盖. for(int k=0;k<n;k++) { Point p[2]; int cnt = LineToCir(l,cs[k],p);//直线与圆的交点. //这一步有问题。 if(cnt==2)//直线必须与圆有两个交点。 { cnt=0; int flag=0; Point tp[5]; if( LEQ( Dis(l.p1,cs[k].c),cs[k].r ) )//p1在圆内 { tp[flag++]=l.p1; } if( LEQ( Dis(l.p2,cs[k].c),cs[k].r ) )//p2在圆内 { tp[flag++]=l.p2; } if( flag != 2 ) { //有问题。 if( p[0]!=l.p1 && p[0]!=l.p2 && OnSeg(p[0],l) ) tp[flag++]=p[0]; if( p[1]!=l.p1 && p[1]!=l.p2 && OnSeg(p[1],l) ) tp[flag++]=p[1]; } if(flag==2) { double dis1,dis2; dis1 = Dis(tp[0],l.p1); dis2 = Dis(tp[1],l.p1); if(dis1>dis2) swap(dis1,dis2); ls[ lcnt ].l=dis1; ls[ lcnt ].r=dis2; lcnt++; } } } // 得出线段与所有圆,所交线段。 sort(ls,ls+lcnt,cmp); double pre=0; int sign=0; for(int k=0;k<lcnt;k++) { if( GT( ls[k].l,pre ) ) { sign=1; break; } pre = max(pre , ls[k].r); } if(sign==0 && EQ(pre,Dis(l.p1,l.p2)) )//EQ必须要 { mat[i][j] = mat[j][i] = Dis(l.p1,l.p2); } } } //接下来就是求一次最短路了。 double ans=Dij(mat,0,1,id); printf("Case %d: ",tt++); if( EQ(ans,-1) ) printf("No such path.\n"); else printf("%.4lf\n",FIX(ans)+EPS); } return 0; }
标签:
原文地址:http://www.cnblogs.com/chenhuan001/p/5121987.html