Interstellar … Fantasy 三维计算几何，简单性质

Interstellar … Fantasy 三维计算几何，简单性质

题意

给定空间中的两点\(A,B\)，以及一个球\((O,R)\)，保证两点不在球内部，求点\(A\)到\(B\)且不经过球的最短路径长

分析

空间中的性质不太好分析，我们知道球心和\(A,B\)确定了一个平面，因此可以转化到平面上考虑

若线段不经过球，显然答案就是两点距离

否则，显然最短长度就是对圆做切线，切线长加圆弧，这部分计算是简单的

判断线段是否与球交，可以判断球心到线段的距离和球半径

于是套上计算几何板子即可，注意板子需要特判两点是同个点的情况

代码

int sgn(double x){
    if(fabs(x) < eps) return 0;
    if(x < 0)return -1;
    return 1;
}

struct Point{
    double x,y,z;
    Point(){}
    Point(double _x,double _y,double _z):x(_x),y(_y),z(_z){}
    void input(){
        x = rd();
        y = rd();
        z = rd();
    }
    double len(){
        return sqrt(x * x + y * y + z * z);
    }
    Point operator - (const Point&b)const{
        return Point(x - b.x,y - b.y,z - b.z);
    }
    double operator * (const Point&b)const{
        return x * b.x + y * b.y + z * b.z;
    }
    Point operator ^(const Point&b)const{
        return Point(y * b.z - z * b.y,z * b.x - x * b.z,x * b.y - y * b.x);
    }
};


inline double dis(Point a,Point b){
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y) + (a.z - b.z) * (a.z - b.z));
}

inline double DIS(Point a,Point b){
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y) + (a.z - b.z) * (a.z - b.z);
}

struct Line{
    Point s,e;
    Line(){}
    Line(Point _s,Point _e):s(_s),e(_e){}
    double length(){
        return dis(s,e);
    }
    double disPtoL(Point p){
        return ((e - s) ^ (p - s)).len() / dis(s,e);
    }
    double disPtoS(Point p){
        if(sgn((p - s) * (e - s)) < 0 || sgn((p - e) * (s - e)) < 0) 
            return min(dis(p,s),dis(p,e));
        return disPtoL(p);
    }
};

int main(){
    int T = rd();
    while(T--){
        Point o;
        o.input();
        double r = rd();
        Point s,t;
        s.input();
        t.input();
        Line l(s,t);
        double D = l.disPtoS(o);
        if(equals(s.x,t.x) && equals(s.y,t.y) && equals(s.z,t.z)) {
            printf("0.000000000\n");
            continue;
        }
        if(D >= r) {
            printf("%.12Lf\n",dis(s,t));
        }
        else {
            double x = asin(r / dis(o,s));
            double l = sqrt(DIS(s,o) - r * r); 
            double xx = asin(r / dis(o,t));
            double ll = sqrt(DIS(t,o) - r * r);
            double d = dis(s,t);
            double a = dis(o,s);
            double b = dis(o,t);
            double theta = acos((a * a + b * b - d * d) / (2 * a * b));
            theta -= acos(r / dis(o,s)) + acos(r / dis(o,t));
            printf("%.12Lf\n",l + ll + r * theta);
        }
    }
}