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);
        }
    }
}