csu 1503: 点到圆弧的距离
1503: 点到圆弧的距离
Time Limit: 1 Sec Memory Limit: 128 MB Special JudgeSubmit: 614 Solved: 101
[Submit][Status][Web Board]
Description
输入一个点P和一条圆弧（圆周的一部分），你的任务是计算P到圆弧的最短距离。换句话说，你需要在圆弧上找一个点，到P点的距离最小。
提示：请尽量使用精确算法。相比之下，近似算法更难通过本题的数据。
Input
输入包含最多10000组数据。每组数据包含 8个整数x1, y1, x2, y2, x3, y3, xp, yp。圆弧的起点是A(x1,y1)，经过点B(x2,y2)，结束位置是C(x3,y3)。点P的位置是 (xp,yp)。输入保证A, B, C各不相同且不会共线。上述所有点的坐标绝对值不超过20。
Output
对于每组数据，输出测试点编号和P到圆弧的距离，保留三位小数。你的输出和标准输出之间最多能有0.001的误差。
Sample Input
0 0 1 1 2 0 1 -1
3 4 0 5 -3 4 0 1
Sample Output
Case 1: 1.414
Case 2: 4.000
HINT
Source
CSU后台出了问题,找了标准输入输出文件,找不出差异。测试数据下载
这个题磨死我了....学长们比赛的时候是怎么AC的...
说下我的思路:我的思路是判断 弧ABC 为劣弧还是优弧,怎么判断 弧ABC是劣弧还是优弧呢?我们可以判断扇形 SAOC 是否等于 SBOC+SAOB 如果相等,那么 ABC 就是劣弧,否则就是优弧.这里要特判一下平角,不然结果会有问题。如果 弧ABC 为劣弧,我们只要求出射线OP和圆的交点P1,那么只要 SAOP1+SP1OC = SAOC，那么P就在弧ABC里面,最短距离就是|OP-r|.否则,P就没有在弧ABC内,最短距离就是min(AP,CP).如果弧ABC为优弧,同样求出SAOC，SAOP1+SP1OC如果 SAOP1+SP1OC = SAOC，那么P就不在弧ABC里面,最短距离就是min(AP,CP).否则,P就在弧ABC内,最短距离就是|OP-r|.如果是角AOC是平角,那么我们就要特判.判断OP和OB是否是否在OC的同一侧(这里可以用叉积判断),在同一侧，距离就是|OP-r|,否则就是min(AP,CP)...还有这个题最坑的一点,幸亏有测试用例,当我们求圆心角时,我是用 a = acos(OA*OB / |OA*OB|) ，但是由于精度原因,acos()有可能不属于 [-1,1]，但是即使误差小,但是计算机也是会算错的,所以就算是acos(1.00000000000000000200)这样求出的结果也是NAN。所以要强制的转换成1，-1.这样的话整个程序就完成了.
#include <cstdio> #include <cstring> #include <algorithm> #include <cmath> #include <map> #define eps (1e-4) #define N 205 #define dd double #define sqr(x) ((x)*(x)) const double pi = acos(-1); using namespace std; struct Point { double x,y; }; double cross(Point a,Point b,Point c) ///叉积 { return (a.x-c.x)*(b.y-c.y)-(b.x-c.x)*(a.y-c.y); } double dis(Point a,Point b) ///距离 { return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); } double mult(Point a,Point b,Point c) ///点积 { return (a.x-c.x)*(b.x-c.x)+(a.y-c.y)*(b.y-c.y); } Point waixin(Point a,Point b,Point c) ///外接圆圆心坐标 { Point p; double a1 = b.x - a.x, b1 = b.y - a.y, c1 = (a1*a1 + b1*b1)/2; double a2 = c.x - a.x, b2 = c.y - a.y, c2 = (a2*a2 + b2*b2)/2; double d = a1*b2 - a2*b1; p.x = a.x + (c1*b2 - c2*b1)/d, p.y=a.y + (a1*c2 -a2*c1)/d; return p; } Point intersection(Point a,Point b,Point c,Point d) { Point p = a; double t = ((a.x-c.x)*(c.y-d.y)-(a.y-c.y)*(c.x-d.x))/((a.x-b.x)*(c.y-d.y)-(a.y-b.y)*(c.x-d .x)); p.x +=(b.x-a.x)*t; p.y +=(b.y-a.y)*t; return p; } void intersection_line_circle(Point c,double r,Point l1,Point l2,Point& p1,Point& p2) { Point p=c; double t; p.x+=l1.y-l2.y; p.y+=l2.x-l1.x; p=intersection(p,c,l1,l2); t=sqrt(r*r-dis(p,c)*dis(p,c))/dis(l1,l2); p1.x=p.x+(l2.x-l1.x)*t; p1.y=p.y+(l2.y-l1.y)*t; p2.x=p.x-(l2.x-l1.x)*t; p2.y=p.y-(l2.y-l1.y)*t; } double getarea(Point a,Point b,Point c,double r){ ///求扇形面积 double temp = mult(a,b,c)/(dis(a,c)*dis(b,c)); if(temp<=-1) temp = -1; ///这里一定要记得判范围..不然 1.00000000000000000200 得到的acos()也是 NAN if(temp>=1) temp = 1; double angle = acos(temp); return fabs(angle*r*r/2); } int main() { Point a,b,c,p,p1,p2,p3; int t = 1; //freopen("a.in","r",stdin); //freopen("a.txt","w",stdout); while(scanf("%lf%lf%lf%lf%lf%lf%lf%lf",&a.x,&a.y,&b.x,&b.y,&c.x,&c.y,&p.x,&p.y)!=EOF) { Point circle = waixin(a,b,c); double r = dis(circle,a); double ans = min(dis(p,a),dis(p,c)); double op = dis(circle,p); double area1 = getarea(a,c,circle,r); double area2 = getarea(a,b,circle,r); double area3 = getarea(b,c,circle,r); intersection_line_circle(circle,r,circle,p,p1,p2); if(acos(mult(p,p1,circle)/((dis(p,circle))*dis(p1,circle)))<0) p3 = p2; else p3 = p1; double area4 = getarea(a,p3,circle,r); double area5 = getarea(c,p3,circle,r); double temp = mult(a,c,circle)/(dis(a,circle)*dis(c,circle)); if(temp<=-1) temp = -1; if(temp>=1) temp = 1; double angle_aoc = acos(temp); if(fabs(angle_aoc-pi)<eps){ ///平角特殊处理 if(cross(p,c,circle)*cross(b,c,circle)>=0){ ans = min(ans,fabs(op-r)); } printf("Case %d: %.3lf\n",t++,ans); continue; } if(fabs(area1-area2-area3)<eps) ///劣弧 { if(fabs(area1-area4-area5)<eps) ans = min(ans,fabs(op-r)); } else { if(fabs(area1-area4-area5)>eps) ans = min(ans,fabs(op-r)); } printf("Case %d: %.3lf\n",t++,ans); } }
解法二:学校队友提供(利用弧度求解)
#include <iostream> #include <cstdio> #include <cstring> #include <cmath> #include <algorithm> #define PI acos(-1.0) using namespace std; double r,ans,dis; double ant,ant1,ant2,ant3; struct Point { double x,y; }p[20]; double jisuan(Point a,Point b) { return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); } bool is_ni() { if(ant3>ant1) { if(ant2<ant3&&ant2>ant1)return false; else return true; } if(ant3<ant1) { if(ant2<ant1&&ant2>ant3)return true; return false; } } void get_fang() { if(p[3].y>0) { if(p[3].x==0)ant=PI/2; else ant=atan2(p[3].y,p[3].x); } else if(p[3].y==0) { if(p[3].x>0)ant=0; else ant=PI; } else if(p[3].y<0) { if(p[3].x==0)ant=3*PI/2; ant=2*PI+atan2(p[3].y,p[3].x); } if(p[0].y>0) { if(p[0].x==0)ant1=PI/2; else ant1=atan2(p[0].y,p[0].x); } else if(p[0].y==0) { if(p[0].x>0)ant1=0; else ant1=PI; } else if(p[0].y<0) { if(p[0].x==0)ant1=3*PI/2; ant1=2*PI+atan2(p[0].y,p[0].x); } if(p[2].y>0) { if(p[2].x==0)ant2=PI/2; else ant2=atan2(p[2].y,p[2].x); } else if(p[2].y==0) { if(p[2].x>0)ant2=0; else ant2=PI; } else if(p[2].y<0) { if(p[2].x==0)ant2=3*PI/2; ant2=2*PI+atan2(p[2].y,p[2].x); } if(p[1].y>0) { if(p[1].x==0)ant3=PI/2; else ant3=atan2(p[1].y,p[1].x); } else if(p[1].y==0) { if(p[1].x>0)ant3=0; else ant3=PI; } else if(p[1].y<0) { if(p[1].x==0)ant3=3*PI/2; ant3=2*PI+atan2(p[1].y,p[1].x); } } Point getyuanxin(Point p1, Point p2, Point p3) { double a, b, c, d, e, f; Point p; a = 2*(p2.x-p1.x); b = 2*(p2.y-p1.y); c = p2.x*p2.x+p2.y*p2.y-p1.x*p1.x-p1.y*p1.y; d = 2*(p3.x-p2.x); e = 2*(p3.y-p2.y); f = p3.x*p3.x+p3.y*p3.y-p2.x*p2.x-p2.y*p2.y; p.x = (b*f-e*c)/(b*d-e*a); p.y = (d*c-a*f)/(b*d-e*a); r = sqrt((p.x-p1.x)*(p.x-p1.x)+(p.y-p1.y)*(p.y-p1.y)); return p; } int main() { // freopen("a.in","r",stdin); // freopen("out.txt","w",stdout); int l=1; while(~scanf("%lf%lf%lf%lf%lf%lf%lf%lf",&p[0].x,&p[0].y,&p[1].x,&p[1].y,&p[2].x,&p[2].y,&p[3].x,&p[3].y)) { int flag=0,flag1=0; p[4]=getyuanxin(p[0],p[1],p[2]); dis=jisuan(p[4],p[3]); r=jisuan(p[4],p[0]); for(int i=0;i<=4;i++)p[i].x-=p[4].x,p[i].y-=p[4].y; get_fang(); if(is_ni()) { if(ant1>ant2) { if(ant>=ant1||ant<=ant2)ans=abs(dis-r); else ans=min(jisuan(p[0],p[3]),jisuan(p[2],p[3])); } else { if(ant<=ant2&&ant>=ant1)ans=abs(dis-r); else ans=min(jisuan(p[0],p[3]),jisuan(p[2],p[3])); } } else { if(ant1>ant2) { if(ant>=ant2&&ant<=ant1)ans=abs(dis-r); else ans=min(jisuan(p[0],p[3]),jisuan(p[2],p[3])); } else { if(ant<=ant1||ant>=ant2)ans=abs(dis-r); else ans=min(jisuan(p[0],p[3]),jisuan(p[2],p[3])); } } printf("Case %d: %.3lf\n",l++,ans); } return 0; }
解法三:经过冥思苦想,终于想到了怎样判断一个点是否在两个向量的夹角之间!
大家请看上图,假设当前 OA 在 OC 的逆时针方向,现在B点绕C点旋转相对于O点是顺时针方向,B点绕A点的话相对于O点是逆时针方向,然后这样就可以判断B点是否在OA~OB之间了，然后马上就可以判断当前弧是优弧还是劣弧,这样的话P点也以同样的方式判断就可以得到其位于优弧还是劣弧了(只不过P点要判0)，但是还是要特判平角.不过平角的判断的判断我也做了一个小小的优化,不用角度判了,只用 2*r == dis(a,c) 判断就OK，这样就不怕被角度坑了.
#include <cstdio> #include <cstring> #include <algorithm> #include <cmath> #include <map> #define eps (1e-4) #define N 205 #define dd double #define sqr(x) ((x)*(x)) const double pi = acos(-1); using namespace std; struct Point { double x,y; }; double cross(Point a,Point b,Point c) ///叉积 { return (a.x-c.x)*(b.y-c.y)-(b.x-c.x)*(a.y-c.y); } double dis(Point a,Point b) ///距离 { return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); } Point waixin(Point a,Point b,Point c) ///外接圆圆心坐标 { Point p; double a1 = b.x - a.x, b1 = b.y - a.y, c1 = (a1*a1 + b1*b1)/2; double a2 = c.x - a.x, b2 = c.y - a.y, c2 = (a2*a2 + b2*b2)/2; double d = a1*b2 - a2*b1; p.x = a.x + (c1*b2 - c2*b1)/d, p.y=a.y + (a1*c2 -a2*c1)/d; return p; } bool judge(Point a,Point b,Point c,Point p){ ///此模板判断点 p 是否在两条射线 ab 和 ac 之间(但是p点不在 ab或者 ac上) if(cross(b,c,a)>0&&cross(p,a,c)<0&&cross(p,a,b)>0){ ///b 在 c 的逆时针方向 return true; } if(cross(b,c,a)<0&&cross(p,a,c)>0&&cross(p,a,b)<0){ ///b 在 c 的顺时针方向 return true; } return false; } bool judge1(Point a,Point b,Point c,Point p){ ///此模板判断点 p 是否在两条射线 ab 和 ac 之间(p点可以在 ab或者 ac上) if(cross(b,c,a)>0&&cross(p,a,c)<=0&&cross(p,a,b)>=0){ ///b 在 c 的逆时针方向 return true; } if(cross(b,c,a)<0&&cross(p,a,c)>=0&&cross(p,a,b)<=0){ ///b 在 c 的顺时针方向 return true; } return false; } int main() { Point a,b,c,p; int t = 1; freopen("a.in","r",stdin); freopen("a.txt","w",stdout); while(scanf("%lf%lf%lf%lf%lf%lf%lf%lf",&a.x,&a.y,&b.x,&b.y,&c.x,&c.y,&p.x,&p.y)!=EOF) { Point circle = waixin(a,b,c); double r = dis(circle,a); double ans = min(dis(p,a),dis(p,c)); double op = dis(circle,p); if(fabs(2*r-dis(a,c))<eps){ ///平角特殊处理 if(cross(p,c,circle)*cross(b,c,circle)>=0){ ans = min(ans,fabs(op-r)); } printf("Case %d: %.3lf\n",t++,ans); continue; } if(judge(circle,a,c,b)) ///劣弧 { if(judge1(circle,a,c,p)) ans = min(ans,fabs(op-r)); } else { if(!judge(circle,a,c,p)) ans = min(ans,fabs(op-r)); } printf("Case %d: %.3lf\n",t++,ans); } }
官方标程:
// Rujia Liu #include<cmath> #include<cstdio> #include<iostream> using namespace std; const double PI = acos(-1.0); const double TWO_PI = 2 * PI; const double eps = 1e-6; inline double NormalizeAngle(double rad, double center = PI) { return rad - TWO_PI * floor((rad + PI - center) / TWO_PI); } inline int dcmp(double x) { if(fabs(x) < eps) return 0; else return x < 0 ? -1 : 1; } struct Point { double x, y; Point(double x=0, double y=0):x(x),y(y) { } }; typedef Point Vector; inline Vector operator + (Vector A, Vector B) { return Vector(A.x+B.x, A.y+B.y); } inline Vector operator - (Point A, Point B) { return Vector(A.x-B.x, A.y-B.y); } inline double Dot(Vector A, Vector B) { return A.x*B.x + A.y*B.y; } inline double Cross(Vector A, Vector B) { return A.x*B.y - A.y*B.x; } inline double Length(Vector A) { return sqrt(Dot(A, A)); } // 外接圆圆心。假定三点不共线 Point get_circumscribed_center(Point p1, Point p2, Point p3) { double bx = p2.x - p1.x; double by = p2.y - p1.y; double cx = p3.x - p1.x; double cy = p3.y - p1.y; double d = 2 * (bx * cy - by * cx); Point p; p.x = (cy * (bx * bx + by * by) - by * (cx * cx + cy * cy)) / d + p1.x; p.y = (bx * (cx * cx + cy * cy) - cx * (bx * bx + by * by)) / d + p1.y; return p; } double DistanceToArc(Point a, Point start, Point mid, Point end) { ///点到圆弧的距离 Point p = get_circumscribed_center(start, mid, end); bool CCW = dcmp(Cross(mid - start, end - start)) > 0; double ang_start = atan2(start.y-p.y, start.x-p.x); double ang_end = atan2(end.y-p.y, end.x-p.x); double r = Length(p - start); double ang = atan2(a.y-p.y, a.x-p.x); bool inside; if(CCW) { inside = NormalizeAngle(ang - ang_start) < NormalizeAngle(ang_end - ang_start); } else { inside = NormalizeAngle(ang - ang_end) < NormalizeAngle(ang_start - ang_end); } if(inside) { return fabs(r - Length(p - a)); } return min(Length(a - start), Length(a - end)); } int main() { int kase = 0; double x1, y1, x2, y2, x3, y3, xp, yp; while(cin >> x1 >> y1 >> x2 >> y2 >> x3 >> y3 >> xp >> yp) { double ans = DistanceToArc(Point(xp,yp), Point(x1,y1), Point(x2,y2), Point(x3,y3)); printf("Case %d: %.3lf\n", ++kase, ans); } return 0; }
