POJ 2177 Ghost Busters(三维几何)
Description
![](http://poj.org/images/2177_1.jpg)
Ghosts are detected by ectoplasmic scanner and are represented as floating spheres. The coordinates of their centers and radii are delivered from the ectoplasmic scanner to the targeting software. The coordinate system is aligned is such a way, that the proton gun fires from the point (0, 0, 0) anywhere into X ≥ 0, Y ≥ 0, Z ≥ 0 trihedral angle. The gun fires a proton ray in a straight line and is so powerful, that even a touch of its ray is enough to kill a ghost. The ray of the proton gun is able to kill a virtually unlimited number of ghosts on its way.
For the first prototype for the targeting software, you are asked to write a program that determines the maximal number of ghosts that can be killed with a single shot of the proton gun.
Input
On the first line of the input there is a single integer number N (0 ≤ N ≤ 100) - the number of ghosts detected by the ectoplasmic scanner of Ecto-1. The following N lines describe detected ghosts - one ghost per line. The description of ith ghost (ghosts are numbered from 1 to N) consists of 4 integer numbers Xi, Yi, Zi, and Ri, separated by spaces. Xi, Yi, Zi (1 ≤ Xi, Yi, Zi ≤ 10000) are the coordinates of the ghost's center, and Ri (1 ≤ Ri ≤ min(Xi, Yi, Zi)) is the ghost's radius. Because ghosts are ectoplasmic, they can be arbitrarily placed in respect to each others. They can intersect, fit inside each other, coincide with each other, etc.
Output
On the first line of the output write a single integer number - the maximal number of ghosts that can be killed with a single shot of the proton gun. On the second line of the output file write the identifying numbers of the ghosts to be killed in an arbitrary order separated by spaces. If there are multiple ways to kill this number of ghosts then write any one.
题目大意:在(x≥0,y≥0,z≥0)的空间直角坐标系上有N个球,问从原点发出一条射线,问最多可以穿过多少个球(接触球的边缘也算),输出任意一组穿过的球。
思路:从原点出发的射线要能穿过两个球,那么,我们从原点看过去,这两个球应该是相交的,但是实际上却不一定相交。
那么,我们可以把球缩放到离原点相同的距离上,令看上去相交的球实际上也相交。
然后,枚举所有射线可能的方向,判断这个射线穿过了多少个球。
枚举的方向分别为:所有球的球心、从原点看过去的这些球的交点(射线的方向应该可以选择一个多个圆的交,那么它的每个顶点应该是由圆的交点所组成的)。
代码(172MS):
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include <cstdio> 2 #include <algorithm> 3 #include <cstring> 4 #include <iostream> 5 #include <cmath> 6 using namespace std; 7 typedef long long LL; 8 9 const double EPS = 1e-6; 10 const double INF = 1e50; 11 const double PI = acos(-1.0); 12 13 inline int sgn(double x) { 14 return (x > EPS) - (x < -EPS); 15 } 16 17 inline double zero(double x) { 18 if(sgn(x) == 0) return 0; 19 else return x; 20 } 21 22 inline double sqr(double x) { 23 return x * x; 24 } 25 26 struct Point3D { 27 double x, y, z; 28 Point3D() {} 29 Point3D(double x, double y, double z): x(x), y(y), z(z) {} 30 void read() { 31 scanf("%lf%lf%lf", &x, &y, &z); 32 } 33 double operator * (const Point3D &rhs) const { 34 return x * rhs.x + y * rhs.y + z * rhs.z; 35 } 36 Point3D operator + (const Point3D &rhs) const { 37 return Point3D(x + rhs.x, y + rhs.y, z + rhs.z); 38 } 39 Point3D operator - (const Point3D &rhs) const { 40 return Point3D(x - rhs.x, y - rhs.y, z - rhs.z); 41 } 42 Point3D operator * (double rhs) const { 43 return Point3D(x * rhs, y * rhs, z * rhs); 44 } 45 Point3D operator / (double rhs) const { 46 return Point3D(x / rhs, y / rhs, z / rhs); 47 } 48 bool operator == (const Point3D &rhs) const { 49 return sgn(x - rhs.x) == 0 && sgn(y - rhs.y) == 0 && sgn(z - rhs.z) == 0; 50 } 51 double length() const { 52 return sqrt(x * x + y * y + z * z); 53 } 54 Point3D unit() const { 55 return *this / length(); 56 } 57 }; 58 59 struct Line3D { 60 Point3D st, ed; 61 Line3D() {} 62 Line3D(Point3D st, Point3D ed): st(st), ed(ed) {} 63 }; 64 65 struct Plane3D { 66 Point3D a, b, c; 67 Plane3D() {} 68 Plane3D(Point3D a, Point3D b, Point3D c): a(a), b(b), c(c) {} 69 void read() { 70 a.read(), b.read(), c.read(); 71 } 72 }; 73 74 struct Circle3D { 75 Point3D c; 76 double r; 77 Circle3D() {} 78 Circle3D(Point3D c, double r): c(c), r(r) {} 79 void read() { 80 c.read(); 81 scanf("%lf", &r); 82 } 83 }; 84 85 double dist(const Point3D &a, const Point3D &b) { 86 return (a - b).length(); 87 } 88 //叉积 89 Point3D cross(const Point3D &u, const Point3D &v) { 90 Point3D ret; 91 ret.x = u.y * v.z - u.z * v.y; 92 ret.y = u.z * v.x - u.x * v.z; 93 ret.z = u.x * v.y - u.y * v.x; 94 return ret; 95 } 96 //点到直线距离 97 double point_to_line(const Point3D &p, const Line3D &l) { 98 return cross(p - l.st, l.ed - l.st).length() / dist(l.ed, l.st); 99 } 100 //求两直线间的距离 101 double line_to_line(const Line3D u, const Line3D v) { 102 Point3D n = cross(u.ed - u.st, v.ed - v.st); 103 return fabs((u.st - v.st) * n) / n.length(); 104 } 105 //取平面法向量 106 Point3D vector_of_plane(const Plane3D &s) { 107 return cross(s.a - s.b, s.b - s.c); 108 } 109 //判断两直线是否平行 110 bool isParallel(const Line3D &u, const Line3D &v) { 111 return sgn(cross(u.ed - u.st, v.ed - v.st).length()) <= 0; 112 } 113 //判断直线是否与球相交 114 bool isIntersect(const Line3D &l, const Circle3D &cir) { 115 return sgn(point_to_line(cir.c, l) - cir.r) <= 0; 116 } 117 //直线与平面的交点 118 Point3D intersect(const Line3D &l, const Plane3D &s) { 119 Point3D ret = vector_of_plane(s); 120 double t = (ret * (s.a - l.st)) / (ret * (l.ed - l.st)); 121 return l.st + (l.ed - l.st) * t; 122 } 123 //在原点上看,两个球的交点 124 int intersect(const Circle3D &u, const Circle3D &v, Point3D &p1, Point3D &p2) { 125 double d = dist(u.c, v.c); 126 if(u.c == v.c || sgn(d - u.r - v.r) > 0 || sgn(fabs(u.r - v.r) - d) > 0) return 0; 127 double t = (sqr(d) + sqr(u.r) - sqr(v.r)) / (2 * d); 128 Point3D mid = u.c + (v.c - u.c).unit() * t; 129 Point3D vec = cross(mid, v.c - u.c).unit() * sqrt(zero(sqr(u.r) - sqr(t))); 130 p1 = mid + vec; 131 p2 = mid - vec; 132 return 1 + sgn(vec.length()); 133 } 134 135 const int MAXN = 110; 136 137 Circle3D cir[MAXN]; 138 Point3D p[MAXN * MAXN], ansVec; 139 int maxAns, pcnt; 140 int n; 141 142 int count(const Point3D &vec) { 143 int ret = 0; 144 for(int i = 0; i < n; ++i) 145 ret += (sgn(point_to_line(cir[i].c, Line3D(Point3D(0, 0, 0), vec)) - cir[i].r) <= 0); 146 return ret; 147 } 148 149 void output(const Point3D &vec) { 150 bool flag = false; 151 for(int i = 0; i < n; ++i) { 152 if(sgn(point_to_line(cir[i].c, Line3D(Point3D(0, 0, 0), vec)) - cir[i].r) <= 0) { 153 if(flag) putchar(' '); 154 flag = true; 155 printf("%d", i + 1); 156 } 157 } 158 printf("\n"); 159 } 160 161 int main() { 162 scanf("%d", &n); 163 for(int i = 0; i < n; ++i) cir[i].read(); 164 for(int i = 0; i < n; ++i) { 165 double t = 20000 / cir[i].c.length(); 166 cir[i].c = cir[i].c * t; 167 cir[i].r = cir[i].r * t; 168 } 169 pcnt = 0; 170 for(int i = 0; i < n; ++i) 171 for(int j = i + 1; j < n; ++j) pcnt += intersect(cir[i], cir[j], p[pcnt], p[pcnt + 1]); 172 maxAns = 0; 173 for(int i = 0; i < n; ++i) { 174 int t = count(cir[i].c); 175 if(t > maxAns) { 176 maxAns = t; 177 ansVec = cir[i].c; 178 } 179 } 180 for(int i = 0; i < pcnt; ++i) { 181 int t = count(p[i]); 182 if(t > maxAns) { 183 maxAns = t; 184 ansVec = p[i]; 185 } 186 } 187 printf("%d\n", maxAns); 188 output(ansVec); 189 }