POJ 2177 Ghost Busters（三维几何）

Description

The famous Ghost Busters team has decided to upgrade their Ectomobile (aka Ecto-1) with a powerful proton gun and an advanced targeting system. Egon has designed and built all the hardware which consists of ectoplasmic scanner and a proton gun that has two degrees of freedom and can automatically rotate and fire in a 90 degrees trihedral angle. You have been hired to write a prototype for the targeting software. 

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）：

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <iostream>
#include <cmath>
using namespace std;
typedef long long LL;

const double EPS = 1e-6;
const double INF = 1e50;
const double PI = acos(-1.0);

inline int sgn(double x) {
    return (x > EPS) - (x < -EPS);
}

inline double zero(double x) {
    if(sgn(x) == 0) return 0;
    else return x;
}

inline double sqr(double x) {
    return x * x;
}

struct Point3D {
    double x, y, z;
    Point3D() {}
    Point3D(double x, double y, double z): x(x), y(y), z(z) {}
    void read() {
        scanf("%lf%lf%lf", &x, &y, &z);
    }
    double operator * (const Point3D &rhs) const {
        return x * rhs.x + y * rhs.y + z * rhs.z;
    }
    Point3D operator + (const Point3D &rhs) const {
        return Point3D(x + rhs.x, y + rhs.y, z + rhs.z);
    }
    Point3D operator - (const Point3D &rhs) const {
        return Point3D(x - rhs.x, y - rhs.y, z - rhs.z);
    }
    Point3D operator * (double rhs) const {
        return Point3D(x * rhs, y * rhs, z * rhs);
    }
    Point3D operator / (double rhs) const {
        return Point3D(x / rhs, y / rhs, z / rhs);
    }
    bool operator == (const Point3D &rhs) const {
        return sgn(x - rhs.x) == 0 && sgn(y - rhs.y) == 0 && sgn(z - rhs.z) == 0;
    }
    double length() const {
        return sqrt(x * x + y * y + z * z);
    }
    Point3D unit() const {
        return *this / length();
    }
};

struct Line3D {
    Point3D st, ed;
    Line3D() {}
    Line3D(Point3D st, Point3D ed): st(st), ed(ed) {}
};

struct Plane3D {
    Point3D a, b, c;
    Plane3D() {}
    Plane3D(Point3D a, Point3D b, Point3D c): a(a), b(b), c(c) {}
    void read() {
        a.read(), b.read(), c.read();
    }
};

struct Circle3D {
    Point3D c;
    double r;
    Circle3D() {}
    Circle3D(Point3D c, double r): c(c), r(r) {}
    void read() {
        c.read();
        scanf("%lf", &r);
    }
};

double dist(const Point3D &a, const Point3D &b) {
    return (a - b).length();
}
//叉积
Point3D cross(const Point3D &u, const Point3D &v) {
    Point3D ret;
    ret.x = u.y * v.z - u.z * v.y;
    ret.y = u.z * v.x - u.x * v.z;
    ret.z = u.x * v.y - u.y * v.x;
    return ret;
}
//点到直线距离
double point_to_line(const Point3D &p, const Line3D &l) {
    return cross(p - l.st, l.ed - l.st).length() / dist(l.ed, l.st);
}
//求两直线间的距离
double line_to_line(const Line3D u, const Line3D v) {
    Point3D n = cross(u.ed - u.st, v.ed - v.st);
    return fabs((u.st - v.st) * n) / n.length();
}
//取平面法向量
Point3D vector_of_plane(const Plane3D &s) {
    return cross(s.a - s.b, s.b - s.c);
}
//判断两直线是否平行
bool isParallel(const Line3D &u, const Line3D &v) {
    return sgn(cross(u.ed - u.st, v.ed - v.st).length()) <= 0;
}
//判断直线是否与球相交
bool isIntersect(const Line3D &l, const Circle3D &cir) {
    return sgn(point_to_line(cir.c, l) - cir.r) <= 0;
}
//直线与平面的交点
Point3D intersect(const Line3D &l, const Plane3D &s) {
    Point3D ret = vector_of_plane(s);
    double t = (ret * (s.a - l.st)) / (ret * (l.ed - l.st));
    return l.st + (l.ed - l.st) * t;
}
//在原点上看，两个球的交点
int intersect(const Circle3D &u, const Circle3D &v, Point3D &p1, Point3D &p2) {
    double d = dist(u.c, v.c);
    if(u.c == v.c || sgn(d - u.r - v.r) > 0 || sgn(fabs(u.r - v.r) - d) > 0) return 0;
    double t = (sqr(d) + sqr(u.r) - sqr(v.r)) / (2 * d);
    Point3D mid = u.c + (v.c - u.c).unit() * t;
    Point3D vec = cross(mid, v.c - u.c).unit() * sqrt(zero(sqr(u.r) - sqr(t)));
    p1 = mid + vec;
    p2 = mid - vec;
    return 1 + sgn(vec.length());
}

const int MAXN = 110;

Circle3D cir[MAXN];
Point3D p[MAXN * MAXN], ansVec;
int maxAns, pcnt;
int n;

int count(const Point3D &vec) {
    int ret = 0;
    for(int i = 0; i < n; ++i)
        ret += (sgn(point_to_line(cir[i].c, Line3D(Point3D(0, 0, 0), vec)) - cir[i].r) <= 0);
    return ret;
}

void output(const Point3D &vec) {
    bool flag = false;
    for(int i = 0; i < n; ++i) {
        if(sgn(point_to_line(cir[i].c, Line3D(Point3D(0, 0, 0), vec)) - cir[i].r) <= 0) {
            if(flag) putchar(' ');
            flag = true;
            printf("%d", i + 1);
        }
    }
    printf("
");
}

int main() {
    scanf("%d", &n);
    for(int i = 0; i < n; ++i) cir[i].read();
    for(int i = 0; i < n; ++i) {
        double t = 20000 / cir[i].c.length();
        cir[i].c = cir[i].c * t;
        cir[i].r = cir[i].r * t;
    }
    pcnt = 0;
    for(int i = 0; i < n; ++i)
        for(int j = i + 1; j < n; ++j) pcnt += intersect(cir[i], cir[j], p[pcnt], p[pcnt + 1]);
    maxAns = 0;
    for(int i = 0; i < n; ++i) {
        int t = count(cir[i].c);
        if(t > maxAns) {
            maxAns = t;
            ansVec = cir[i].c;
        }
    }
    for(int i = 0; i < pcnt; ++i) {
        int t = count(p[i]);
        if(t > maxAns) {
            maxAns = t;
            ansVec = p[i];
        }
    }
    printf("%d
", maxAns);
    output(ansVec);
}