SGU 110 Dungeon(立體幾何)

Description:

The mission of space explorers found on planet M the vast dungeon. One of the dungeon halls is fill with the bright spheres. The explorers find out that the light rays reflect from the surface of the spheres according the ordinary law (the incidence angle is equal to the reflectance angle, the incidence ray, the reflected ray and the perpendicular to the sphere surface lay in the one plane). The ancient legend says that if the light ray will reflect from the spheres in the proper order, than the door to the room with very precious ancient knowledge will open. You are not to guess the right sequence; your task is much simpler. You are given the positions and the radii of the spheres, the place where the laser shot was made and the direction of light propagation. And you must find out the sequence in which the light will be reflected from the spheres.

Input:

The first line of input contains the single integer n (1≤n≤50) - the amount of the spheres. The next n lines contain the coordinates and the radii of the spheres xi, yi, zi, ri (the integer numbers less or equal to 10000 by absolute value). The last line contains 6 real numbers - the coordinates of two points. The first one gives the coordinates of the place of laser shot, and the second gives the direction in which it was made (the second point is the point on the ray). The starting point of the ray lies strictly outside of any sphere.

Output:

Your program must output the sequence of sphere numbers (spheres are numbers from 1 as they was given in input), from which the light ray was reflected. If the ray will reflect more the 10 times, than you must output first 10, then a space and the word ‘etc.’ (without quotes). Notice: if the light ray goes at a tangent to the sphere you must assume that the ray was reflected by the sphere.

Sample Input:

1
0 0 2 1
0 0 0 0 0 1

Sample Output:

1

Sample Input:

2
0 0 2 1
0 0 -2 1
0 0 0 0 0 100

Sample Output:

1 2 1 2 1 2 1 2 1 2 etc.

題目連結

在空間內有 $n$ 個球和一條射線，射線照射到球面上會發生反射（相切也算反射），求射線路徑中依次發生反射的球體編號（超過10次反射僅輸出10次反射的球體編號和etc.）

從題目資料範圍看對於每次的射線可以暴力列舉每個球判斷其是否相交，其中相交切交點最近的球為反射球，再求出反射射線重複上述過程即可

首先需要判斷射線是否與球體相交，顯然射線經球反射後的反射線一定在射線的起點、秋心和射線與球的交點三點確定的平面上，這樣就可以把立體幾何先轉化到平面幾何上然後再繼續分析

若射線為圖中 $CD$ ，球為 $\odot A$ ，射線與球的交點為 $D、F$ ，反射線為 $DC'$ （ $C'$ 為 $C$ 關於法線 $AD$ 的對稱點 ）

我們只討論射線與球的一個交點 $D$ ，$E$ 暫且不討論

根據向量的關係可得 $\vec{DA}=\vec{CA}-\vec{CD}$ 而 $|\vec{DA}|=r$ （ $r$ 是球 $A$ 的半徑）

假設 $|\vec{CD}|=l$ , $\vec{CD}$ 的方向向量為 $\vec{dir}$， $CD=l\times\vec{dir}$

聯立得 $|\vec{CA}-l\times\vec{dir}|=R$ 其中 $l$ 為變數

這樣我們就得到了一個關於 $l$ 的一元二次方程，方程化簡到最後為
$\vec{AC}^{2}+2l\vec{AC}\vec{dir}+l^{2}\vec{dir}^{2}=r^{2}$
移項得
$l^{2}\vec{dir}^{2}+2l\vec{AC}\vec{dir}+\vec{AC}^{2}r^{2}=0$
設 $a=\vec{dir}^{2},b={2\vec{AC}\vec{dir}},c=\vec{AC}^{2}$ ，根據根的判別式 $\Delta=b^{2}-4ac$ 就可以判斷射線所在的直線是否與球體相交，若相交則繼續利用求根公式求出兩解 $\frac{-b\pm\sqrt{\Delta}}{2a}$ ，取最小解（顯然圖中 $|\vec{CD}|<|\vec{CE}|$ ），若最小解為負數則射線所在直線與球體相交但射線不與球體相交（球體在射線反向），若最小解為 $0$ 則射線與球體相切，若最小解為整數則射線與球體相交

判斷完是否相交就需要求出射線經球體反射後的反射線了

因為此題的射線是用兩點（射線起點與射線上另一點）表示，所以反射線的起點可以用入射線與球體的交點表示，而另一反射線上的點利用入射線起點對於法線的對稱點表示就可以啦

AC程式碼:

#include <bits/stdc++.h>
using namespace std;
typedef double db;
const int maxn = 1e2 + 5;
const db inf = 1e20;
const db eps = 1e-9;

int Sgn(db Key) {return fabs(Key) < eps ? 0 : (Key < 0 ? -1 : 1);}
int Cmp(db Key1, db Key2) {return Sgn(Key1 - Key2);}
struct Point {db X, Y, Z;};
typedef Point Vector;
struct Ray {Point Origin; Vector Dir;};
Vector operator - (Vector Key1, Vector Key2) {return (Vector){Key1.X - Key2.X, Key1.Y - Key2.Y, Key1.Z - Key2.Z};}
Vector operator + (Vector Key1, Vector Key2) {return (Vector){Key1.X + Key2.X, Key1.Y + Key2.Y, Key1.Z + Key2.Z};}
db operator * (Vector Key1, Vector Key2) {return Key1.X * Key2.X + Key1.Y * Key2.Y + Key1.Z * Key2.Z;}
db GetLen(Vector Key) {return sqrt(Key * Key);}
db GetLen2(Vector Key) {return Key * Key;}
db operator ^ (Vector Key1, Vector Key2) {return GetLen((Vector){Key1.Y * Key2.Z - Key1.Z * Key2.Y, Key1.Z * Key2.X - Key1.X * Key2.Z, Key1.X * Key2.Y - Key1.Y * Key2.X});}
Vector operator * (Vector Key1, db Key2) {return (Vector){Key1.X * Key2, Key1.Y * Key2, Key1.Z * Key2};}
Vector operator / (Vector Key1, db Key2) {return (Vector){Key1.X / Key2, Key1.Y / Key2, Key1.Z / Key2};}
db DisPointToPoiont(Point Key1, Point Key2) {return GetLen(Key2 - Key1);}
struct Sphere {Point Center; db Radius;};

bool IsRayInterSphere(Ray Key1, Sphere Key2, db &Dis) {
    db A = Key1.Dir * Key1.Dir;
    db B = (Key1.Origin - Key2.Center) * Key1.Dir * 2.0;
    db C = ((Key1.Origin - Key2.Center) * (Key1.Origin - Key2.Center)) - (Key2.Radius * Key2.Radius);
    db Delta = B * B - 4.0 * A * C;
    if (Sgn(Delta) < 0) return false;
    db X1 = (-B - sqrt(Delta)) / (2.0 * A), X2 = (-B + sqrt(Delta)) / (2.0 * A);
    if (Cmp(X1, X2) > 0) swap(X1, X2);
    if (Sgn(X1) <= 0) return false;
    Dis = X1;
    return true;
}

void Reflect(Ray &Key1, Sphere Key2, db Dis) {
    Point Pos = Key1.Origin + (Key1.Dir * Dis);
    Vector Temp = Key2.Center + (((Pos - Key2.Center) * ((Pos - Key2.Center) * (Key1.Origin - Key2.Center))) / GetLen2(Pos - Key2.Center));
    Key1.Dir = Temp * 2.0 - Key1.Origin - Pos; Key1.Origin = Pos;
}

int N;
Sphere spheres[maxn];
Ray Cur;
int Last;

int Solve() {
    db MinDis = inf; int Ans = -1;
    for (int i = 1; i <= N; ++i) {
        db Dis = inf;
        if (i == Last) continue;
        if (!IsRayInterSphere(Cur, spheres[i], Dis)) continue;
        if (Cmp(Dis, MinDis) < 0) {
            Ans = i;
            MinDis = Dis;
        }
    }
    if (Ans == -1) return -1;
    Last = Ans;
    Reflect(Cur, spheres[Ans], MinDis);
    return Ans;
}

int main(int argc, char *argv[]) {
    scanf("%d", &N);
    for (int i = 1; i <= N; ++i) scanf("%lf%lf%lf%lf", &spheres[i].Center.X, &spheres[i].Center.Y, &spheres[i].Center.Z, &spheres[i].Radius);
    scanf("%lf%lf%lf", &Cur.Origin.X, &Cur.Origin.Y, &Cur.Origin.Z);
    scanf("%lf%lf%lf", &Cur.Dir.X, &Cur.Dir.Y, &Cur.Dir.Z);
    Cur.Dir = Cur.Dir - Cur.Origin;
    for (int i = 0, Ans; i < 11; ++i) {
        Ans = Solve();
        if (Ans == -1) break;
        else if (i == 10) printf("etc.");
        else printf("%d ", Ans);
    }
    printf("\n");
    return 0;
}