2009年11月30日星期一.sgu110
sgu110:完全被這題玩了。。。
題意很簡單：
一堆球，一條光線，光線反射符合反射定律，輸出所有反射經(jīng)過的球的序號
的確沒有光線的起點(diǎn)在球中。
tricks:
1. test case 3 是無交點(diǎn)的。。。。我被這個玩了
2. 直線要表示成參數(shù)方程
P = P0 + s*v; 注意是射線，所以只有 "s >= 0" 時是成立的,注意 >= 0,自己想下為什么有等于0
還要特判一下不能連著交同一個圓
貌似沒有別的問題了
反射直線有兩個推法
1.
lt表示直線的參數(shù)方程
lt.p 為p0,lt.v為v
point_t v = normal(tmp - circles[idx]);
double d = -(tmp.x * v.x + tmp.y * v.y + tmp.z * v.z); //平面方程ax + by + cz + d = 0中的d
double dis = fabs(p2p(lt.p, v.x, v.y, v.z, d)); // p2p為點(diǎn)到平面的距離
v = tmp - lt.p + scale(v, dis * 2);
lt.p = tmp, lt.v = normal(v);
下圖表示了上面求解的過程
\ | / \ | / \|/ --------- \ \ \ 還可以跟據(jù),
2.
求出中點(diǎn)，mid = (s1 + s2)/2;
然后求出s2即可
\ * / \ / \ /==============================華麗麗的分割線==============================
貼代碼
1 /*
2 * SOUR:sgu110
3 * ALGO:3D computational geometry
4 * DATE: 2009年 11月 29日 星期日 17:21:22 CST
5 * COMM:4
6 * */
7 #include<iostream>
8 #include<cstdio>
9 #include<cstdlib>
10 #include<cstring>
11 #include<algorithm>
12 #include<cmath>
13 using namespace std;
14 typedef long long LL;
15 const int maxint = 0x7fffffff;
16 const long long max64 = 0x7fffffffffffffffll;
17 const int N = 128;
18 int n;
19 struct point_t {
20 double x,y,z,r;
21 point_t(){}
22 point_t(double a,double b,double c){x = a,y = b,z = c;}
23 }circle[N],s,e;
24 point_t operator + (point_t a,point_t b) { return point_t(a.x + b.x ,a.y + b.y,a.z + b.z);}
25 point_t operator - (point_t a,point_t b) { return point_t(a.x - b.x ,a.y - b.y,a.z - b.z);}
26 point_t scale(point_t a,double fac) { return point_t(a.x * fac,a.y * fac,a.z * fac);}
27 double dist(point_t a) { return sqrt(a.x * a.x + a.y * a.y + a.z * a.z);}
28 point_t normal(point_t a) {
29 double len = dist(a);
30 return point_t(a.x / len,a.y/len,a.z/len);
31 }
32 point_t setLen(point_t a,double fac) {
33 a = normal(a);
34 return point_t(a.x*fac,a.y*fac,a.z*fac);
35 }
36 int out[N*N],top;
37
38 double p2p(point_t p,double a,double b,double c,double d)
39 {
40 return (a * p.x + b * p.y + c * p.z + d)/sqrt(a * a + b * b + c* c);
41 }
42
43 const double inf = 1e10;
44 const double eps = 1e-12;
45 double dot_mul(point_t a,point_t b) {
46 return a.x * b.x + a.y * b.y + a.z * b.z;
47 }
48 void solve()
49 {
50 int i,j,k,pre = -1;
51 top = 0;
52 for(k = 0;k < 12;k++) {
53 int idx = 0;
54 double mink = inf;
55 for(i = 1;i <= n;i++) {
56 if(pre != i) {
57 double a = (e.x * e.x + e.y * e.y + e.z * e.z);
58 point_t t = s - circle[i];
59 double b = 2 * dot_mul(e,t);
60 double c = t.x * t.x + t.y * t.y + t.z * t.z - circle[i].r * circle[i].r;
61 if(b * b >= 4 * a * c) {
62 double s1 = (-b - sqrt(b*b - 4 * a*c))/(2*a);
63 double s2 = (-b + sqrt(b*b - 4 * a*c))/(2*a);
64 if(s1 >= 0) {
65 if(s1 < mink) {
66 mink = s1;
67 idx = i;
68 }
69 }else if(s2 >= 0){
70 if(s2 < mink) {
71 mink = s1;
72 idx = i;
73 }
74 }
75 }
76 }
77 }
78 if(idx > 0) {
79 point_t p = s + scale(e,mink);
80 //p is the point which intersects with the circles
81 point_t v = normal(p - circle[idx]); // normal vector
82 double d = -dot_mul(p,v);
83 //double dis = fabs(p2p(s,v.x,v.y,v.z,d));
84 double dis = fabs(dot_mul(v,s)+d);
85
86 point_t mid = p + scale(v,dis);
87 point_t p1 = scale(mid,2) - s;
88
89 s = p,e =p1-p;
90 out[top++] = idx;
91 pre = idx;
92 }else {
93 break;
94 }
95 }
96 if(top > 0) //哥死在這句話上了。。。
97 printf("%d",out[0]);
98 for(i = 1;i < 10 && i < top;i++) {
99 printf(" %d",out[i]);
100 }
101 if(top > 10) {
102 printf(" etc.");
103 }
104 putchar(10);
105 }
106
107 int main()
108 {
109 int i,j,k;
110 while(scanf("%d",&n) == 1) {
111 top = 0;
112 for(i = 1;i <= n;i++) {
113 scanf("%lf%lf%lf%lf",&circle[i].x,&circle[i].y,&circle[i].z,&circle[i].r);
114 }
115 scanf("%lf%lf%lf%lf%lf%lf",&s.x,&s.y,&s.z,&e.x,&e.y,&e.z);
116 e = e - s;
117 solve();
118 }
119 return 0;
120 }
121