POJ 3944 - Spherical Mirrors

http://poj.org/problem?id=3944

概要

三次元空間に N 個の球がある．この球は図のようにレーザーを反射する．

(0,0,0) から (u,v,w) の方向にレーザーを発射したとき，最後に反射する球面上の点を答える．

制約

解法

三次元幾何をがんばる． ライブラリを用意できていれば特に難しい点は無い．

半直線と球の交差判定は，半直線と円の交差判定と同様に，向きを確認してから点と直線の距離が r 未満かどうかを見ればいい．

半直線と球の交点は二分探索で求めた．

反射させるときは，2回の外積で QN に関して QI とは逆方向の法線を出して，QI と QN の内積で得られた cosθ を使って QR を作った． ただし，θ = 0 のときは QI と QN の外積が零ベクトルになってしまうので，単に入射光のベクトルを反転させたものを返すようにした．

  1 #include <cstdio>
  2 #include <vector>
  3 #include <map>
  4 #include <cmath>
  5 using namespace std;
  6 static const double EPS = 1e-6;
  7 
  8 struct P
  9 {
 10   double x, y, z;
 11   P() {}
 12   P(double a, double b, double c) : x(a), y(b), z(c) {}
 13   P operator+(const P& that) const { return P(x + that.x, y + that.y, z + that.z); }
 14   P operator-(const P& that) const { return P(x - that.x, y - that.y, z - that.z); }
 15   P& operator*=(double a) { x *= a; y *= a; z *= a; return *this; }
 16   P operator*(double a) const { P p(*this); p *= a; return p; }
 17   P& operator/=(double a) { x /= a; y /= a; z /= a; return *this; }
 18   P operator/(double a) const { P p(*this); p /= a; return p; }
 19 
 20   double abs() const { return sqrt(x*x + y*y + z*z); }
 21   double dot(const P& that) const { return x*that.x + y*that.y + z*that.z; }
 22   P cross(const P& that) const
 23   {
 24     return P(
 25         y*that.z - that.y*z,
 26         z*that.x - that.z*x,
 27         x*that.y - that.x*y);
 28   }
 29 };
 30 
 31 struct line
 32 {
 33   P a, b;
 34   line(const P& x, const P& y) : a(x), b(y) {}
 35 };
 36 
 37 struct sphere
 38 {
 39   P o;
 40   double r;
 41   sphere(const P& p, double x) : o(p), r(x) {}
 42 
 43   bool intersects(const line& ln) const
 44   {
 45     const P ab = ln.b - ln.a;
 46     const P ao = o - ln.a;
 47     if (ab.dot(ao) < 0) {
 48       return false;
 49     }
 50     const double h = ab.cross(ao).abs() / ab.abs();
 51     return h < r;
 52   }
 53 
 54   P intersection(const line& ln) const
 55   {
 56     const P ab = ln.b - ln.a;
 57     const P ao = o - ln.a;
 58     const double h = ab.cross(ao).abs() / ab.abs();
 59     const double d = ao.abs();
 60     const double l = sqrt(d*d - h*h);
 61     double left = 0.0, right = l/ab.abs();
 62     for (int i = 0; i < 100; i++) {
 63       const double mid = (left + right)/2.0;
 64       const P ap = ab*mid;
 65       const double t = (ap - ao).abs();
 66       if (t < r) {
 67         right = mid;
 68       } else {
 69         left = mid;
 70       }
 71     }
 72     return ln.a + ab*((left + right)/2.0);
 73   }
 74 
 75   line reflect(const line& ln) const
 76   {
 77     const P q = intersection(ln);
 78     const P qa = ln.a - q;
 79     P qn = q - o;
 80     qn /= qn.abs();
 81     const P x = qn.cross(qa);
 82     if (x.abs() < EPS) {
 83       return line(q, ln.a);
 84     }
 85     P qm = qn.cross(x);
 86     qm /= qm.abs();
 87     const double t = qn.dot(qa)/(qn.abs()*qa.abs());
 88     const P qr = qm*sqrt(1-t*t) + qn*t;
 89     return line(q, q + qr);
 90   }
 91 };
 92 
 93 int main()
 94 {
 95   int N;
 96   while (scanf("%d", &N) != EOF && N != 0) {
 97     int u, v, w;
 98     scanf("%d %d %d", &u, &v, &w);
 99     line ln(P(0, 0, 0), P(u, v, w));
100     vector<sphere> ss;
101     for (int i = 0; i < N; i++) {
102       int x, y, z, r;
103       scanf("%d %d %d %d", &x, &y, &z, &r);
104       ss.push_back(sphere(P(x, y, z), r));
105     }
106 
107     P ans;
108     for (;;) {
109       map<double,int> m;
110       for (int i = 0; i < N; i++) {
111         if (ss[i].intersects(ln)) {
112           const P p = ss[i].intersection(ln);
113           m.insert(make_pair((p - ln.a).abs(), i));
114         }
115       }
116       if (m.empty()) {
117         break;
118       }
119       const sphere& s = ss[m.begin()->second];
120       ln = s.reflect(ln);
121       ans = ln.a;
122     }
123     printf("%.4f %.4f %.4f\n", ans.x, ans.y, ans.z);
124   }
125   return 0;
126 }
poj/3944.cc