hdu4063(圆与圆交+线段与圆交+最短路)
写几何题总是提心吊胆。精度问题真心吓人。
其实思路挺简单的一道题,真是什么算法和几何double搞到一块,心里就虚虚的。
思路:求出所有圆之间的交点,然后用这些交点跑一遍最短路就可以了。
Aircraft
Time Limit: 10000/3000 MS (Java/Others) Memory Limit: 65536/32768 K (Java/Others)
Total Submission(s): 1244 Accepted Submission(s): 304
Problem Description
You are playing a flying game.
In the game, player controls an aircraft in a 2D-space.
The mission is to drive the craft from starting point to terminal point.
The craft needs wireless signal to move.
A number of devices are placed in the 2D-space, spreading signal.
For a device Di, it has a signal radius -- Ri.
When the distance between the craft and Di is shorter or equal to Ri, it(the craft) gets Di's wireless signal.
Now you need to tell me the shortest path from starting point to terminal point.
In the game, player controls an aircraft in a 2D-space.
The mission is to drive the craft from starting point to terminal point.
The craft needs wireless signal to move.
A number of devices are placed in the 2D-space, spreading signal.
For a device Di, it has a signal radius -- Ri.
When the distance between the craft and Di is shorter or equal to Ri, it(the craft) gets Di's wireless signal.
Now you need to tell me the shortest path from starting point to terminal point.
Input
The first line of the input file is a single integer T.
The rest of the test file contains T blocks.
Each block starts with an integer n, followed by n devices given as (xi, yi, Ri).
(xi, yi) is position of Di, and Ri is the radius of its signal range.
The first point is the starting point.
The last point is the terminal point.
T <= 25;
2 <= n <= 20 for most cases;
20 < n <= 25 for several cases, completely random generated.
-1000 <= xi, yi <= 1000 , 1 <= ri <= 1000.
All are integers.
The rest of the test file contains T blocks.
Each block starts with an integer n, followed by n devices given as (xi, yi, Ri).
(xi, yi) is position of Di, and Ri is the radius of its signal range.
The first point is the starting point.
The last point is the terminal point.
T <= 25;
2 <= n <= 20 for most cases;
20 < n <= 25 for several cases, completely random generated.
-1000 <= xi, yi <= 1000 , 1 <= ri <= 1000.
All are integers.
Output
For each case, Output "No such path." if the craft can't get to the terminal point.
Otherwise, output a float number, correct the result to 4 decimal places.(as shown in the sample output)
Otherwise, output a float number, correct the result to 4 decimal places.(as shown in the sample output)
Sample Input
2
2
0 0 1
2 0 1
2
0 0 1
4 1 2
Sample Output
Case 1: 2.0000
Case 2: No such path.
#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; }