UVa11887 Tetrahedrons and Spheres

题目链接

    本题是2010年湖南省第六届大学生计算机程序设计竞赛的K题

题意

    三维空间中有a个四面体和b个球（1≤a,b≤5），求它们体积的并（重叠部分体积只累加一次）。所有坐标为绝对值不超过5的整数，球的半径r_i≤3。

分析

    肯定要用到数值积分，用自适应辛普森法来求，主要还是看按z值切片后的平面问题该怎么处理：如何求三角形、四边形、圆的面积并。​
    如果没有圆，求多边形的面积，可以用梯形剖分的思想来处理：假设多边形的顶点是逆时针顺序给出的，遍历每条边时，如果顶点是从左到右的（横坐标从小到大）就把两个顶点之间与过最低点的水平线组成的一个直角梯形作为负面积，相反从右往左则把构成的直角梯形作为正面积，然后用正面积减掉负面积就可以得到多边形的面积。
    要求所有多边形的面积并，则需要找出所有的梯形分割点：多边形本身的顶点横坐标和多边形相交产生的点的横坐标。分割点排序后，逐段（$[x_{a-1},x_a]$）对多边形切割得出所有的正负梯形，再求面积并。
   细说一下逐段（$[x_{a-1},x_a]$）求面积并的处理过程：先把梯形斜边按位置关系从上到下排序（一般用梯形斜边的中点的高度来比大小），正/负面积梯形有一个符号标识$f_i$（取值$\pm 1$），面积记录一个正值（实际上记录梯形斜边中点的高度$h_i$），初始计数标识$cnt=0$以及面积结果为0，遍历梯形时如果$cnt>0$则面积累加$(x_a-x_{a-1})\times (h_{i-1}-h_i)$，再更新$cnt=cnt+f_i$。
   其实对这种方法稍加改造就能处理有圆的情况：分割点加入圆心以及圆的左右端点、圆与圆交点以及多边形与圆交点的横坐标，再逐段切割多变形和圆，切割圆时同样找出正负梯形部分并求出梯形斜边与圆弧那部分的面积$a_i$，再把梯形斜边中点竖直向外延伸至与圆相交处的高度作为排序依据。同样在遍历时如果$cnt>0$则对面积进行累加，累加量变成：$(x_a-x_{a-1})\times (h_{i-1}-h_i)+a_{i-1}\times f_{i-1}-a_i\times f_i$。
   可以设想一下多边形和圆相交的各种情况，画一个分割后的图模拟并深入理解一下这个过程。
   还有一个更简洁高效也更直观的处理方式：将梯形从上到下排序后，对各连通部分找出最上最下两个梯形就能得出此部分面积为$$(x_a-x_{a-1})\times (h_i-h_j)+a_i+a_j$。
   给一份测试数据。
   最后附两个梯形剖分相关的RMQ问题：
   梯形剖分入门(特殊情况)【ice】
   考试B 冰雪奇缘改版 多边形剖梯形+线段树维护区间

AC 代码

#include <iostream>
#include <iomanip>
#include <cmath>
#include <algorithm>
using namespace std;

#define eps 1e-12
#define M 335
#define N 5
int x[N][4], y[N][4], z[N][4], cx[N], cy[N], cz[N], r[N], c[N], a, b, m, t;
double px[N][5], py[N][5], sx[M], cr[N];

double interp(double x1, double y1, double x2, double y2, double x) {
    return (y1*(x2-x) + y2*(x-x1)) / (x2-x1);
}

void cross(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4) {
    double dx = x2-x1, dy = y2-y1, dx1 = x3-x1, dy1 = y3-y1, dx2 = x4-x1, dy2 = y4-y1;
    if ((dx*dy1-dx1*dy)*(dx*dy2-dx2*dy) >= 0.) return;
    dx = x4-x3; dy = y4-y3; dx1 = x1-x3; dy1 = y1-y3; dx2 = x2-x3; dy2 = y2-y3;
    double c1 = dx*dy1-dx1*dy, c2 = dx*dy2-dx2*dy;
    if (c1*c2 >= 0.) return;
    c1 /= c1 - c2; sx[t++] = (1-c1)*x1 + c1*x2;
}

void cross(double x1, double y1, double x2, double y2, int i) {
    double dx = x2-x1, dy = y2-y1, l = sqrt(dx*dx + dy*dy);
    if (l < eps) return;
    double vx = cx[i]-x1, vy = cy[i]-y1; dx /= l; dy /= l;
    double d = abs(dx*vy-vx*dy);
    if (d > cr[i]-eps) return;
    double p = dx*vx+dy*vy, q = sqrt(cr[i]*cr[i]-d*d), p1 = p-q, p2 = p+q;
    if (p1 > -eps && p1 < l+eps) sx[t++] = x1 + dx*p1;
    if (p2 > -eps && p2 < l+eps) sx[t++] = x1 + dx*p2;
}

void cross(int i, int j) {
    double dx = cx[i]-cx[j], dy = cy[i]-cy[j], d = sqrt(dx*dx + dy*dy);
    if (d > cr[i]+cr[j]-eps || d < abs(cr[i]-cr[j])+eps) return;
    double cosa = dx/d, cosb = (cr[j]*cr[j]+d*d-cr[i]*cr[i])/2./cr[j]/d, ss = dy/d*sqrt(1-cosb*cosb), cc = cosa*cosb;
    sx[t++] = cx[j]+cr[j]*(cc-ss); sx[t++] = cx[j]+cr[j]*(cc+ss);
}

struct trap {
    double h, t, a; int f;
    bool operator< (const trap& rhs) const {
        return t > rhs.t;
    }
} s[N<<2];

void trap(int i, double l, double r) {
    int cc = 0; double h[2];
    for (int j=0; j<c[i]; ++j) {
        double x1 = px[i][j], y1 = py[i][j], x2 = px[i][j+1], y2 = py[i][j+1];
        if (min(x1,x2) > l+eps || max(x1,x2) < r-eps) continue;
        h[cc++] = .5 * (interp(x1, y1, x2, y2, l) + interp(x1, y1, x2, y2, r));
    }
    if (cc != 2) return;
    if (h[0] < h[1]) h[0] += h[1], h[1] = h[0]-h[1], h[0] = h[0]-h[1];
    s[m].h = h[0]; s[m].t = h[0]; s[m].a = 0.; s[m++].f = 1; s[m].h = h[1]; s[m].t = h[1]; s[m].a = 0.; s[m++].f = -1;
}

void trapc(int i, double l, double r) {
    if (max(abs(l-cx[i]), abs(r-cx[i])) > cr[i]+eps) return;
    double d1 = cr[i]*cr[i]-(l-cx[i])*(l-cx[i]), d2 = cr[i]*cr[i]-(r-cx[i])*(r-cx[i]);
    d1 = d1<=0. ? 0. : sqrt(d1); d2 = d2<=0. ? 0. : sqrt(d2);
    double x = .5 * (l+r), d = cr[i]*cr[i]-(x-cx[i])*(x-cx[i]); d = d<=0. ? 0. : sqrt(d);
    double ds = ((l-r)*(l-r) + (d1-d2)*(d1-d2)) / 4., h = sqrt(cr[i]*cr[i]-ds); ds = sqrt(ds);
    double a = cr[i]*cr[i]*asin(ds/cr[i]) - ds*h;
    s[m].h = cy[i] + .5 * (d1 + d2); s[m].t = cy[i] + d; s[m].a = a; s[m++].f = 1;
    s[m].h = cy[i] - .5 * (d1 + d2); s[m].t = cy[i] - d; s[m].a = a; s[m++].f = -1;
}

double area(double h) {
    for (int i=t=0; i<a; ++i) {
        for (int j=c[i]=0; j<4; ++j) {
            if (z[i][j] == h) {
                px[i][c[i]] = x[i][j]; py[i][c[i]++] = y[i][j]; continue;
            }
            for (int k=j+1; k<4; ++k) if (z[i][k] != h) {
                int e = min(z[i][j], z[i][k]), f = max(z[i][j], z[i][k]);
                if (e==f || f<h || e>h) continue;
                px[i][c[i]] = interp(z[i][j], x[i][j], z[i][k], x[i][k], h);
                py[i][c[i]++] = interp(z[i][j], y[i][j], z[i][k], y[i][k], h);
            }
        }
        if (c[i] < 3) continue;
        for (int j=0; j<c[i]; ++j) sx[t++] = px[i][j];
        if (c[i] == 4) {
            double dx = px[i][2] - px[i][0], dy = py[i][2] - py[i][0],
                dx1 = px[i][1] - px[i][0], dy1 = py[i][1] - py[i][0],
                dx2 = px[i][3] - px[i][0], dy2 = py[i][3] - py[i][0], p;
            if ((dx*dy1-dx1*dy)*(dx*dy2-dx2*dy) > 0.) {
                if ((dx1*dy-dx*dy1)*(dx1*dy2-dx2*dy1) > 0.) {
                    p = px[i][1]; px[i][1] = px[i][2]; px[i][2] = px[i][3]; px[i][3] = p;
                    p = py[i][1]; py[i][1] = py[i][2]; py[i][2] = py[i][3]; py[i][3] = p;
                } else
                    p = px[i][2], px[i][2] = px[i][1], px[i][1] = p, p = py[i][2], py[i][2] = py[i][1], py[i][1] = p;
            }
        }
        px[i][c[i]] = px[i][0]; py[i][c[i]] = py[i][0];
    }
    for (int i=0; i<b; ++i) {
        cr[i] = r[i]*r[i] - (cz[i]-h)*(cz[i]-h); cr[i] = cr[i]<=0. ? 0. : sqrt(cr[i]);
        if (cr[i] <= eps) continue;
        sx[t++] = cx[i]; sx[t++] = cx[i]-cr[i]; sx[t++] = cx[i]+cr[i];
    }
    for (int i=0; i<a; ++i) if (c[i] > 2) for (int j=0; j<c[i]; ++j) {
        for (int k=i+1; k<a; ++k) if (c[k] > 2) for (int p=0; p<c[k]; ++p)
            cross(px[i][j], py[i][j], px[i][j+1], py[i][j+1], px[k][p], py[k][p], px[k][p+1], py[k][p+1]);
        for (int k=0; k<b; ++k) if (cr[k] > eps) cross(px[i][j], py[i][j], px[i][j+1], py[i][j+1], k);
    }
    for (int i=0; i<b; ++i) if (cr[i] > eps) for (int j=i+1; j<b; ++j) if (cr[j] > eps) cross(i, j);
    sort(sx, sx+t);
    double ans = 0., d;
    for (int i=1; i<t; ++i) if ((d = sx[i]-sx[i-1]) > eps) {
        for (int j=m=0; j<a; ++j) if (c[j] > 2) trap(j, sx[i-1], sx[i]);
        for (int j=0; j<b; ++j) if (cr[j] > eps) trapc(j, sx[i-1], sx[i]);
        sort(s, s+m);
        for (int j=0, cc=0, p; j<m; ++j) {
            if (cc == 0) p = j;
            if ((cc += s[j].f) == 0) ans += d * (s[p].h-s[j].h) + s[p].a + s[j].a;
        }
    }
    return ans;
}

double asr(double l, double m, double r, double fl, double fm, double fr, double ans, double err, int step) {
    double lm = (l+m)/2., flm = area(lm), sl = (fl + 4.*flm + fm) * (m-l) / 6.,
        rm = (m+r)/2., frm = area(rm), sr = (fm + 4.*frm + fr) * (r-m) / 6.;
    if (abs(sl+sr-ans) <= 15.*err && step < 0) return sl + sr + (sl+sr-ans) / 15.;
    return asr(l, lm, m, fl, flm, fm, sl, err/2., step-1) + asr(m, rm, r, fm, frm, fr, sr, err/2., step-1);
}

double asr(double l, double r, double err) {
    double m = (l+r)/2., fl = area(l), fm = area(m), fr = area(r);
    return asr(l, m, r, fl, fm, fr, (fl + 4.*fm + fr) * (r-l) / 6., err, 6);
}

void solve() {
    for (int i=0; i<a; ++i) for (int j=0; j<4; ++j) cin >> x[i][j] >> y[i][j] >> z[i][j];
    for (int i=0; i<b; ++i) cin >> cx[i] >> cy[i] >> cz[i] >> r[i];
    cout << asr(-8., 8., 1e-3) << endl;
}

int main() {
    cout << fixed << setprecision(3);
    while (cin >> a >> b && (a || b)) solve();
    return 0;
}