POJ 2149 - Inherit the Spheres
http://poj.org/problem?id=2149
概要
三次元空間に N 個の球がある．
ある z0 に対して，z = z0 平面において，いくつの物体があるかに興味がある． z0 を 0 から 36000 に動かしたとき，物体の個数の変化の様子を
- 増えたときに 1
- 減ったときに 0
で表したものを答える．
制約
- 1 <= N <= 100
- 1 <= 半径 R[i] <= 200
- 0 < X[i]-R[i] < X[i]+R[i] < 4000
- 0 < Y[i]-R[i] < Y[i]+R[i] < 16000
- 0 < Z[i]-R[i] < Z[i]+R[i] < 36000
- 球と球が1点で接していることはない
- Z[i]±R[i] および球同士の交差によって生じる円のz座標の最大値・最小値は互いに少なくとも 0.01 異なる
注意
- ある球が他の球を含んでいることがある
解法
個数が変わる可能性があるのは，制約にあるように
- Z[i]±R[i]
- 球同士の交差によって生じる円のz座標の最大値・最小値
の2種類なので，これらの z 座標をソートして各段階で何個の物体があるかどうか数えて変化を調べる．
まず球同士の交線について． これを真面目に方程式で解こうとするとしんどいので，二分探索で求めることにする． 球の中心を結ぶ線分上で，どこが円の中心 z0 になるかを方程式から求めて，この z0 から上方向と下方向にそれぞれ， その z で球を切ったときにちょうど切断面の円と円が接するところを二分探索で求める．
次にある z における物体の数について． z で球を切って生じた円に対し，それぞれ交差判定をして Union-Find で数えればいい．
最初は
- Z[i]+R[i] で -1
- Z[i]-R[i] で +1
- 交線の上側で +1
- 交線の下側で -1
として，同じ z 座標になるようなものをうまくマージしてイベントドリブン風に解く方法を考えていた． しかしこれだと3つ以上の球が交わっているようなとき，実際には +1 も -1 もされないような場合があって，それを弾く方法に自信がなかった．
poj/2149.cc1 #include <iostream> 2 #include <vector> 3 #include <algorithm> 4 #include <complex> 5 #include <cmath> 6 #include <iterator> 7 using namespace std; 8 typedef complex<double> P; 9 static const double EPS = 1e-6; 10 11 struct P3 12 { 13 double x, y, z; 14 P3() {} 15 P3(double a, double b, double c) : x(a), y(b), z(c) {} 16 P3& operator+=(const P3& q) { x += q.x; y += q.y; z += q.z; return *this; } 17 P3& operator-=(const P3& q) { x -= q.x; y -= q.y; z -= q.z; return *this; } 18 P3& operator*=(double a) { x *= a; y *= a; z *= a; return *this; } 19 }; 20 P3 operator+(const P3& p, const P3& q) { P3 r(p); return r += q; } 21 P3 operator-(const P3& p, const P3& q) { P3 r(p); return r -= q; } 22 P3 operator*(const P3& p, double a) { P3 q(p); q *= a; return q; } 23 P3 operator*(double a, const P3& p) { return p*a; } 24 inline double dot(const P3& p, const P3& q) { return p.x*q.x + p.y*q.y + p.z*q.z; } 25 inline double abs(const P3& p) { return sqrt(dot(p, p)); } 26 27 struct circle 28 { 29 P o; 30 double r; 31 circle() {} 32 circle(const P& p, double x) : o(p), r(x) {} 33 34 inline bool contains(const circle& c) const { return abs(o - c.o)+c.r <= r; } 35 inline bool intersects(const circle& c) const { return !contains(c) && !c.contains(*this) && abs(o - c.o) <= r + c.r; } 36 pair<P,P> intersection(const circle& c) const 37 { 38 // assert(intersects(c)) 39 const double d = abs(o - c.o); 40 const double cos_ = (d*d + r*r - c.r*c.r) / (2*d); 41 const double sin_ = sqrt(r*r - cos_*cos_); 42 const P e = (c.o - o) / d; 43 return make_pair(o + e*P(cos_, sin_), o + e*P(cos_, -sin_)); 44 } 45 }; 46 47 struct sphere 48 { 49 P3 o; 50 double r; 51 sphere(const P3& p, double x) : o(p), r(x) {} 52 53 inline circle cutz(double z) const { return circle(P(o.x, o.y), sqrt(r*r - (z-o.z)*(z-o.z))); } 54 inline bool contains(const sphere& s) const { return abs(o - s.o)+s.r <= r; } 55 inline bool intersects(const sphere& s) const { return !contains(s) && !s.contains(*this) && abs(o - s.o) <= r + s.r; } 56 57 double center_z(const sphere& s) const 58 { 59 const P3 p = s.o - o; 60 const double L = abs(p); 61 const double x1 = (r*r - s.r*s.r + L*L)/(2*L); 62 return (o + x1/L*p).z; 63 } 64 65 pair<double,double> intersection_z(const sphere& that) const 66 { 67 const double z0 = center_z(that); 68 double left = z0, right = min(o.z+r, that.o.z+that.r); 69 for (int i = 0; i < 100; i++) { 70 const double mid = (left+right)/2.0; 71 const circle c1 = cutz(mid), c2 = that.cutz(mid); 72 const double d = abs(c1.o - c2.o); 73 if (d < c1.r + c2.r) { 74 left = mid; 75 } else { 76 right = mid; 77 } 78 } 79 const double z1 = (left+right)/2.0; 80 81 left = max(o.z-r, that.o.z-that.r), right = z0; 82 for (int i = 0; i < 100; i++) { 83 const double mid = (left+right)/2.0; 84 const circle c1 = cutz(mid), c2 = that.cutz(mid); 85 const double d = abs(c1.o - c2.o); 86 if (d < c1.r + c2.r) { 87 right = mid; 88 } else { 89 left = mid; 90 } 91 } 92 const double z2 = (left+right)/2.0; 93 return make_pair(z1, z2); 94 } 95 }; 96 97 struct DisjointSet/*{{{*/ 98 { 99 vector<int> parent; 100 101 int root(int x) 102 { 103 if (parent[x] < 0) { 104 return x; 105 } else { 106 parent[x] = root(parent[x]); 107 return parent[x]; 108 } 109 } 110 111 explicit DisjointSet(int n) : parent(n, -1) {} 112 113 bool unite(int x, int y) 114 { 115 const int a = root(x); 116 const int b = root(y); 117 if (a != b) { 118 if (parent[a] < parent[b]) { 119 parent[a] += parent[b]; 120 parent[b] = a; 121 } else { 122 parent[b] += parent[a]; 123 parent[a] = b; 124 } 125 return true; 126 } else { 127 return false; 128 } 129 } 130 131 bool find(int x, int y) { return root(x) == root(y); } 132 int size(int x) { return -parent[root(x)]; } 133 };/*}}}*/ 134 135 int count_sphere(const vector<sphere>& ss, double z) 136 { 137 vector<circle> cs; 138 for (vector<sphere>::const_iterator it = ss.begin(); it != ss.end(); ++it) { 139 if (abs(z - it->o.z) < it->r) { 140 cs.push_back(it->cutz(z)); 141 } 142 } 143 144 const int N = cs.size(); 145 DisjointSet s(N); 146 for (int i = 0; i < N; i++) { 147 for (int j = i+1; j < N; j++) { 148 if (cs[i].contains(cs[j]) || cs[j].contains(cs[i]) || cs[i].intersects(cs[j])) { 149 s.unite(i, j); 150 } 151 } 152 } 153 int ans = 0; 154 for (int i = 0; i < N; i++) { 155 if (s.parent[i] < 0) { 156 ++ans; 157 } 158 } 159 return ans; 160 } 161 162 int main() 163 { 164 int N; 165 while (cin >> N && N != 0) { 166 vector<sphere> ss; 167 vector<double> zs; 168 for (int i = 0; i < N; i++) { 169 int x, y, z, r; 170 cin >> x >> y >> z >> r; 171 ss.push_back(sphere(P3(x, y, z), r)); 172 zs.push_back(z-r); 173 zs.push_back(z+r); 174 } 175 176 for (vector<sphere>::const_iterator it = ss.begin(); it != ss.end(); ++it) { 177 for (vector<sphere>::const_iterator jt = it+1; jt != ss.end(); ++jt) { 178 if (it->intersects(*jt)) { 179 const pair<double,double> r = it->intersection_z(*jt); 180 zs.push_back(r.first); 181 zs.push_back(r.second); 182 } 183 } 184 } 185 sort(zs.begin(), zs.end()); 186 187 vector<char> ans; 188 int prev = 0; 189 for (vector<double>::const_iterator it = zs.begin(); it != zs.end(); ++it) { 190 const int r = count_sphere(ss, *it+EPS); 191 if (prev < r) { 192 ans.push_back('1'); 193 } else if (prev > r) { 194 ans.push_back('0'); 195 } 196 prev = r; 197 } 198 199 cout << ans.size() << endl; 200 copy(ans.begin(), ans.end(), ostream_iterator<char>(cout, "")); 201 cout << endl; 202 } 203 return 0; 204 }