toj 4615 Tetrahedrons and Spheres

时间限制(普通/Java):6000MS/18000MS 内存限制:65536KByte
总提交: 2 测试通过:2

描述

There are a tetrahedrons and b spheres in the 3D-splace, you’re asked to calculate the volume occupied by at least one of them (i.e. volume of the union of the objects).

输入

There will be at most 20 test cases. Each case begins with two integers a, b, the number of tetrahedrons and the number of spheres (1<=a,b<=5). The next a lines each contains 12 integers: x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, the coordinates (xi, yi, zi)(1<=i<=4) of the four vertices of a tetrahedron. The next b lines each contains 4 integers x, y, z, r, the coordinates of the center (x, y, z) and the radius r (r<=3). All the coordinate values are integers with absolute values no more than 5. The input is terminated by a=b=0.

输出

For each test case, print a single line, the volume occupied by at least one of them, rounded to three decimal points.

样例输入

1 1
0 0 4 1 0 4 0 1 4 0 0 5
0 0 0 1
0 0

样例输出

4.356

参考程序

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#include <vector>using namespace std;
#define sqr(x) ((x)*(x))const double eps = 1e-6;
const double inf = 8;const double pi = 3.14159265358979323846;inline bool cmpDouble(const double &a, const double &b) { return fabs(a - b) < eps; }struct Tpoint {double x, y;Tpoint() {}Tpoint(double a, double b) {x = a;y = b;}inline double norm() { return sqrt(sqr(x) + sqr(y)); }
};inline Tpoint operator+(const Tpoint &a, const Tpoint &b) { return Tpoint(a.x + b.x, a.y + b.y); }inline Tpoint operator-(const Tpoint &a, const Tpoint &b) { return Tpoint(a.x - b.x, a.y - b.y); }inline Tpoint operator*(const double &a, const Tpoint &b) { return Tpoint(a * b.x, a * b.y); }inline Tpoint operator*(const Tpoint &b, const double &a) { return Tpoint(a * b.x, a * b.y); }inline Tpoint operator/(const Tpoint &a, const double &b) { return Tpoint(a.x / b, a.y / b); }inline bool operator<(const Tpoint &a, const Tpoint &b) {return a.x + eps < b.x || fabs(a.x - b.x) < eps && a.y + eps < b.y;
}inline bool operator==(const Tpoint &a, const Tpoint &b) { return fabs(a.x - b.x) < eps && fabs(a.y - b.y) < eps; }inline double det(const Tpoint &a, const Tpoint &b) { return a.x * b.y - a.y * b.x; }inline double dot(const Tpoint &a, const Tpoint &b) { return a.x * b.x + a.y * b.y; }struct P3 {double x, y, z;P3() {}P3(double a, double b, double c) {x = a;y = b;z = c;}inline void read() { scanf("%lf%lf%lf", &x, &y, &z); }
};inline P3 operator+(const P3 &a, const P3 &b) { return P3(a.x + b.x, a.y + b.y, a.z + b.z); }inline P3 operator-(const P3 &a, const P3 &b) { return P3(a.x - b.x, a.y - b.y, a.z - b.z); }inline P3 operator*(const double &a, const P3 &b) { return P3(a * b.x, a * b.y, a * b.z); }inline P3 operator*(const P3 &b, const double &a) { return P3(a * b.x, a * b.y, a * b.z); }inline P3 operator/(const P3 &a, const double &b) { return P3(a.x / b, a.y / b, a.z / b); }struct Tcir {double r;Tpoint o;Tcir() {}Tcir(Tpoint x, double y) { o = x, r = y; }
};vector<Tcir> Circle;
typedef vector<Tpoint> TP;
vector<TP> Poly;struct Tinter {double x, y, Area, mid;int delta;Tinter() {}Tinter(double xx, double yy, double mm, int dd, double aa) {x = xx;y = yy;mid = mm;delta = dd;Area = aa;}
};inline bool operator<(const Tinter &a, const Tinter &b) { return a.mid > b.mid + eps; }inline bool operator==(const Tinter &a, const Tinter &b) { return fabs(a.mid - b.mid) < eps; }vector<Tinter> inter;
vector<double> bak;inline double dist(const Tpoint &a, const Tpoint &b) {return sqr(a.x - b.x) + sqr(a.y - b.y);
}inline void Add(double x) {bak.push_back(x);
}inline void CircleIntersectCircle(const Tcir &a, const Tcir &b) {double l = dist(a.o, b.o);double s = ((a.r - b.r) * (a.r + b.r) / l + 1) / 2;double t = sqrt(-(l - sqr(a.r + b.r)) * (l - sqr(a.r - b.r)) / (l * l * 4));double ux = b.o.x - a.o.x, uy = b.o.y - a.o.y;double ix = a.o.x + s * ux + t * uy, iy = a.o.y + s * uy - t * ux;double jx = a.o.x + s * ux - t * uy, jy = a.o.y + s * uy + t * ux;Add(ix);Add(jx);
}inline bool InLine(const Tpoint &a, const Tpoint &b, const Tpoint &c) {return fabs(det(b - a, c - a)) < eps && dot(a - c, b - c) < eps;
}inline void LineToLine(const Tpoint &a, const Tpoint &b, const Tpoint &c, const Tpoint &d) {double s1 = det(c - a, b - a), s2 = det(b - a, d - a);if (s1 * s2 < -eps) return;Tpoint e = c + (d - c) * s1 / (s1 + s2);if (InLine(a, b, e) && InLine(c, d, e)) {Add(e.x);}
}inline void LineToCircle(const Tpoint &a, const Tpoint &b, const Tcir &c) {double h = fabs(det(c.o - a, b - a)) / (b - a).norm();if (h > c.r + eps) return;double lamda = dot(c.o - a, b - a);lamda /= dist(a, b);Tpoint x = a + (b - a) * lamda;double d = sqrt(sqr(c.r) - sqr(h));d /= (b - a).norm();Tpoint e = x + (b - a) * d;Tpoint f = x - (b - a) * d;if (InLine(a, b, e))Add(e.x);if (InLine(a, b, f))Add(f.x);return;
}inline void CircleToCircle() {for (int i = 0; i < Circle.size(); ++i) {Add(Circle[i].o.x - Circle[i].r);Add(Circle[i].o.x + Circle[i].r);Add(Circle[i].o.x);for (int j = i + 1; j < Circle.size(); ++j)if (dist(Circle[i].o, Circle[j].o) <= sqr(Circle[i].r + Circle[j].r)) if (dist(Circle[i].o, Circle[j].o) >=sqr(Circle[i].r - Circle[j].r))CircleIntersectCircle(Circle[i], Circle[j]);}
}inline void CircleToPoly() {for (int i = 0; i < Circle.size(); ++i)for (int j = 0; j < Poly.size(); ++j)for (int v = 0; v < Poly[j].size(); ++v)LineToCircle(Poly[j][v], Poly[j][(v + 1) % Poly[j].size()], Circle[i]);
}inline void PolyToPoly() {for (int i = 0; i < Poly.size(); ++i)for (int u = 0; u < Poly[i].size(); ++u)Add(Poly[i][u].x);for (int i = 0; i < Poly.size(); ++i)for (int j = i + 1; j < Poly.size(); ++j)for (int u = 0; u < Poly[i].size(); ++u)for (int v = 0; v < Poly[j].size(); ++v)LineToLine(Poly[i][u], Poly[i][(u + 1) % Poly[i].size()], Poly[j][v],Poly[j][(v + 1) % Poly[j].size()]);
}inline void Get(const Tcir &c, double x, double &l, double &r) {double y = fabs(c.o.x - x);double d = sqrt(fabs(sqr(c.r) - sqr(y)));l = c.o.y + d;r = c.o.y - d;
}inline double arcArea(const Tcir &a, double l, double x, double r, double y) {double len = sqrt(sqr(l - r) + sqr(x - y));double d = sqrt(sqr(a.r) - sqr(len) / 4.0);double angle = atan(len / 2.0 / d);return fabs(angle * sqr(a.r) - d * len / 2.0);
}inline void Get_Interval(const Tcir &a, double l, double r) {double L1, L2, R1, R2, M1, M2;Get(a, l, L1, L2);Get(a, r, R1, R2);Get(a, (l + r) / 2.0, M1, M2);int D1 = 1, D2 = -1;double A1 = arcArea(a, l, L1, r, R1), A2 = arcArea(a, l, L2, r, R2);inter.push_back(Tinter(L1, R1, M1, D1, A1));inter.push_back(Tinter(L2, R2, M2, D2, A2));
}inline double calcSlice(double xl, double xr) {inter.clear();double lmost = -inf, rmost = inf;for (int i = 0; i < Poly.size(); ++i) {int cc = 0;Tinter I[5];for (int u = 0; u < Poly[i].size(); ++u) {Tpoint x = Poly[i][u];Tpoint y = Poly[i][(u + 1) % Poly[i].size()];double l = min(x.x, y.x), r = max(x.x, y.x);if (l <= xl + eps && xr <= r + eps) {if (fabs(l - r) < eps) continue;Tpoint d = y - x;Tpoint Left = x + d / d.x * (xl - x.x);Tpoint Right = x + d / d.x * (xr - x.x);Tpoint Mid = (Left + Right) / 2;I[cc++] = Tinter(Left.y, Right.y, Mid.y, 1, 0);}}sort(I, I + cc);if (cc == 2) {I[1].delta = -1;inter.push_back(I[0]);inter.push_back(I[1]);lmost = max(lmost, I[1].mid);rmost = min(rmost, I[0].mid);}}for (int i = 0; i < Circle.size(); ++i)if (fabs(Circle[i].o.x - xl) < Circle[i].r + eps && fabs(Circle[i].o.x - xr) < Circle[i].r + eps)Get_Interval(Circle[i], xl, xr);if (!inter.size()) return 0;double ans = 0;sort(inter.begin(), inter.end());int cnt = 0;for (int i = 0; i < inter.size(); ++i) {if (cnt > 0) {ans += (fabs(inter[i - 1].x - inter[i].x) + fabs(inter[i - 1].y - inter[i].y)) * (xr - xl) / 2.0;ans += inter[i - 1].delta * inter[i - 1].Area;ans -= inter[i].delta * inter[i].Area;}cnt += inter[i].delta;}return ans;
}#define maxn 105
int n, m;
struct Poly4 {P3 p[4];
} a[maxn];struct Sphere {P3 o;double r;inline void read() {o.read();scanf("%lf", &r);}
} b[maxn];inline bool equal(double a, double b) {return fabs(a - b) < eps;
}inline bool InterSect(const Tpoint &a, const Tpoint &b, const Tpoint &c, const Tpoint &d) {double s1 = det(b - a, c - a), s2 = det(b - a, d - a);if (s1 * s2 > -eps) return false;s1 = det(d - c, a - c), s2 = det(d - c, b - c);if (s1 * s2 > -eps) return false;return true;
}inline void ToHull(vector<Tpoint> &a) {sort(a.begin(), a.end());int hull[10], len, limit = 1;hull[len = 1] = 0;for (int i = 1; i < 4; ++i) {while (len > limit && det(a[hull[len]] - a[hull[len - 1]], a[i] - a[hull[len]]) >= 0) --len;hull[++len] = i;}limit = len;for (int i = 2; i >= 0; --i) {while (len > limit && det(a[hull[len]] - a[hull[len - 1]], a[i] - a[hull[len]]) >= 0) --len;hull[++len] = i;}vector<Tpoint> b = a;a.resize(len - 1);for (int i = 0; i < len - 1; ++i)a[i] = b[hull[i + 1]];
}inline double calcArea(double z) {Poly.clear();Circle.clear();bak.clear();for (int i = 0; i < n; ++i) {vector<Tpoint> cross;for (int j = 0; j < 4; ++j)for (int k = j + 1; k < 4; ++k)if (!equal(a[i].p[j].z, a[i].p[k].z)) {double l = min(a[i].p[j].z, a[i].p[k].z), r = max(a[i].p[j].z, a[i].p[k].z);if (l <= z + eps && z <= r + eps) {P3 d = a[i].p[k] - a[i].p[j];d = d / d.z;d = d * (z - a[i].p[j].z);d = d + a[i].p[j];Tpoint x(d.x, d.y);cross.push_back(x);}}sort(cross.begin(), cross.end());cross.erase(unique(cross.begin(), cross.end()), cross.end());if (cross.size() > 2) {if (cross.size() > 4) {puts("Too Many Points!!!");while (1);}if (cross.size() == 4)ToHull(cross);Poly.push_back(cross);}}for (int i = 0; i < m; ++i)if (fabs(z - b[i].o.z) + eps < b[i].r) {Tpoint o(b[i].o.x, b[i].o.y);double r = sqrt(sqr(b[i].r) - sqr(z - b[i].o.z));Circle.push_back(Tcir(o, r));}CircleToCircle();CircleToPoly();PolyToPoly();sort(bak.begin(), bak.end());double res = 0;for (int i = 0; i + 1 < bak.size(); ++i)if (fabs(bak[i + 1] - bak[i]) > eps)res += calcSlice(bak[i], bak[i + 1]);return res;
}int main() {for (; scanf("%d%d", &n, &m) != EOF && (n + m);) {for (int i = 0; i < n; ++i)for (int j = 0; j < 4; ++j)a[i].p[j].read();for (int i = 0; i < m; ++i)b[i].read();double last = 0, Ans = calcArea(-inf) + calcArea(inf);const int Block = 4000;double h = (inf + inf) / (double) Block;for (int i = 1; i < Block; i += 2)Ans += 4 * calcArea(-inf + i * h);for (int i = 2; i < Block; i += 2)Ans += 2 * calcArea(-inf + i * h);Ans = Ans * h / 3.0;printf("%.3f\n", Ans);}return 0;
}