• csu 1503: 点到圆弧的距离


    1503: 点到圆弧的距离

    Time Limit: 1 Sec  Memory Limit: 128 MB  Special Judge
    Submit: 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
    ",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
    ",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
    ",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
    ",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
    ",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
    ", ++kase, ans);
      }
      return 0;
    }
  • 相关阅读:
    【Thinkphp教程】URL路由功能解析
    MYSQL 错误#145解决方法
    【Thinkphp教程】空模块
    【Thinkphp教程】 如何进行模块分组
    mySQL中删除unique key的语法
    使用php让浏览器刷新
    Spring+Jpa整合的过程中遇到的一个问题。。。纠结了我半天。。。
    关于mule studio的应用
    解决eclipse和myeclipse不能编译项目的问题
    ajax fileupload上传组件的使用感悟
  • 原文地址:https://www.cnblogs.com/liyinggang/p/5775505.html
Copyright © 2020-2023  润新知