UVa1065/LA3809 Raising the Roof

题目链接

  本题是2007年icpc世界总决赛的H题

题意

  给定一些三角形的屋顶，你的任务是计算出暴露在阳光下的部分（即从天空向地面俯视时能看到的部分）的总面积。注意，被遮挡的部分不应该被算到总面积中。xy平面是地面，z轴竖直向上。

输入格式

  输入包含多组数据。输入的第一行为顶点总数V和三角形总数T（1≤V≤300，1≤T≤1000）；以下V行每行3个整数x、y 和z，即各个顶点的坐标；再以下T行每行3个整数，即3个顶点的编号（各顶点按照输入顺序编号为1~V）。所有坐标均为不超过100的正整数。输入没有退化的三角形。任意两个输入三角形的内部没有公共点（边界可以接触，但不会部分重合或者相交）。输入结束标志为V=0且T=0。

输出格式

  对于每组数据，输出可见部分的总面积，保留两位小数。

分析

  和求圆的面积并的经典题目UVa1367/LA3532 Nuclear Plants一个套路，用离散扫描法来做。详细思路《算法竞赛入门经典–训练指南》第435页有写，不再展开叙述。
  本题在vjudge尝试的人数很少，通过的更少，原始oj（UVa）上尝试并通过的人也不多。早期的World Final题目并不给出测试数据，特地开一篇博文填补Raising the Roof测试数据和可视化的坑。

测试数据生成

  按题面交代：没有退化的三角形。任意两个输入三角形的内部没有公共点（边界可以接触，但不会部分重合或者相交）。生成测试数据要满足这些限制，需要一些三维向量运算、平面和三维直线相关的知识储备。

  说一些三维直线相关的知识：空间直线是两个平面相交而来，因此其一般式方程为$\{\begin{array}{l}{a}_{1}x+b_{1}y+c_{1}z+d1=0\\ a_{2}x+b_{2}y+c_{2}z+d2=0\end{array}$。两平面的法向量分别为$\overrightarrow{v_1}=(a_1,b_1,c_1),\overrightarrow{v_2}=(a_2,b_2,c_2)$，因此直线的方向向量$\vec{d}=\overrightarrow{v_1}\times \overrightarrow{v_2}$，如果再求出直线上的任意一点，则可以得到直线的点向式方程。

  计算得到方向向量数值为$\vec{d}=(a,b,c)$，可以这样选取及计算直线上的任意一点$p_0=(x_0,y_0,z_0)$的坐标：
  若$a\ne 0$，则可取$x_0=1$，代入其一般式方程$\{\begin{array}{l}{a}_1+b_1{y}_0+c_1{z}_0+d1=0\\ a_2+b_2{y}_0+c_2{z}_0+d2=0\end{array}$可以计算出$y_0,z_0$；
  若$b\ne 0$，则可取$y_0=1$，代入一般式方程可以计算出$x_0,z_0$；
  若$c\ne 0$，则可取$z_0=1$，代入一般式方程可以计算出$x_0,y_0$。

  其他三维向量运算及立体几何的知识不展开细说，给出一份完整的生成测试数据的python脚本：

# -*- coding: utf-8 -*-from random import sample, randint
from fractions import Fraction as fracX, V, T, C = 100, 10, 15, 10
# X, V, T, C = 100, 5, 7, 300def gen_point(p):while True:x, y, z, ok = randint(1, X), randint(1, X), randint(1, X), Truefor xi, yi, zi in p:if xi == x and yi == y and zi == z:ok = Falsebreakif ok:return x, y, zdef cross(v1, v2):return (v1[1]*v2[2] - v1[2]*v2[1], v1[2]*v2[0] - v1[0]*v2[2], v1[0]*v2[1] - v1[1]*v2[0])def dot(v1, v2):return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]def sub(v1, v2):return (v1[0] - v2[0], v1[1] - v2[1], v1[2] - v2[2])def solve_equations(a1, b1, c1, a2, b2, c2):'''a1*x + b1*y = c1''''''a2*x + b2*y = c2'''p = a1*b2 - a2*b1return frac(b2*c1 - b1*c2, p), frac(a1*c2 - a2*c1, p)def line_intersect_seg(p, n, a, b):v1, v2 = sub(p, a), sub(b, a)x, y, z = cross(n, v2)if x == 0 and y == 0 and z == 0:x, y, z = cross(v1, v2)return [a, b] if x == 0 and y == 0 and z == 0 else []s = (x**2 + y**2 + z**2)**.5x, y, z = cross(n, v1)s = (x**2 + y**2 + z**2)**.5 / sif s >= 0. and s <= 1.: ### eps?return [(a[0] + v2[0]*s, a[1] + v2[1]*s, a[2] + v2[2]*s)]else:return []def line_intersect_tri(p, n, a, b, c):q = line_intersect_seg(p, n, a, b)for x in line_intersect_seg(p, n, a, c):q.append(x)for x in line_intersect_seg(p, n, b, c):q.append(x)if len(q) == 0:return []q.sort(key = lambda x: dot(sub(x, p), n))x0, x1 = q[0], q[-1]return [x0] if x0[0] == x1[0] and x0[1] == x1[1] and x0[2] == x1[2] else [x0, x1]def line_intvs_intersect(a, b, intv):v, d = sub(b, a), []s = v[0]**2 + v[1]**2 + v[2]**2for x in intv:d.append(dot(sub(x, a), v) / s)d.sort()s0, s1 = max(0., d[0]), min(1., d[-1])if s0 > s1:return []return [(a[0] + v[0]*s0, a[1] + v[1]*s0, a[2] + v[2]*s0), (a[0] + v[0]*s1, a[1] + v[1]*s1, a[2] + v[2]*s1)]def point_in_tri(p, a, b, c):c1, c2, c3 = cross(sub(b, a), sub(p, a)), cross(sub(c, b), sub(p, b)), cross(sub(a, c), sub(p, c))return dot(c1, c2) > 0 and dot(c2, c3) > 0 and dot(c3, c1) > 0def seg_cross_seg(a, b, c, d, ec = False):v1, v2, v3 = sub(b, a), sub(d, c), sub(c, a)c1, c2, c3 = cross(v1, v2), cross(v3, v2), cross(v3, v1)if c1[0] == 0 and c1[1] == 0 and c1[2] == 0:return Falsel1 = (c1[0]**2 + c1[1]**2 + c1[2]**2)**.5l2, l3 = (c2[0]**2 + c2[1]**2 + c2[2]**2)**.5 / l1, (c3[0]**2 + c3[1]**2 + c3[2]**2)**.5 / l1if l2 < 0. or l2 > 1.: ### eps?return Falsereturn l3 >= 0. and l3 <= 1. if ec else l3 > 0. and l3 < 1.def seg_cross_tri(a, b, c, d, e):if seg_cross_seg(a, b, c, d):return seg_cross_seg(a, b, c, e, True) or seg_cross_seg(a, b, d, e, True)if seg_cross_seg(a, b, c, e):return seg_cross_seg(a, b, c, d, True) or seg_cross_seg(a, b, d, e, True)if seg_cross_seg(a, b, d, e):return seg_cross_seg(a, b, c, d, True) or seg_cross_seg(a, b, c, e, True)return Falsedef valid(t1, t2):v1, v2, v3, v4 = sub(t1[1], t1[0]), sub(t1[2], t1[0]), sub(t2[1], t2[0]), sub(t2[2], t2[0])n1, n2 = cross(v1, v2), cross(v3, v4)if n1[0] == 0 and n1[1] == 0 and n1[2] == 0:return Falsecc, d1 = cross(n1, n2), dot(n1, t1[0])if cc[0] == 0 and cc[1] == 0 and cc[2] == 0:if dot(n1, t2[0]) != d1:return Truefor p in t1:if point_in_tri(p, t2[0], t2[1], t2[2]):return Falsefor p in t2:if point_in_tri(p, t1[0], t1[1], t1[2]):return Falsefor i in range(1, 3):for j in range(i):if seg_cross_tri(t1[i], t1[j], t2[0], t2[1], t2[2]) or seg_cross_tri(t2[i], t2[j], t1[0], t1[1], t1[2]):return Falsereturn Trueelse:x, y, z, d2 = 1, 1, 1, dot(n2, t2[0])if cc[0] != 0:y, z = solve_equations(n1[1], n1[2], d1 - n1[0]*x, n2[1], n2[2], d2 - n2[0]*x)elif cc[1] != 0:x, z = solve_equations(n1[0], n1[2], d1 - n1[1]*y, n2[0], n2[2], d2 - n2[1]*y)else:x, y = solve_equations(n1[0], n1[1], d1 - n1[2]*z, n2[0], n2[1], d2 - n2[2]*z)q1, q2 = line_intersect_tri((x, y, z), cc, t1[0], t1[1], t1[2]), line_intersect_tri((x, y, z), cc, t2[0], t2[1], t2[2])if len(q1) * len(q2) < 2:return Trueq = line_intvs_intersect(q1[0], q1[1], q2) if len(q1) > 1 else line_intvs_intersect(q2[0], q2[1], q1)if len(q) == 0:return Truep = ((q[0][0]+q[-1][0]) / 2, (q[0][1]+q[-1][1]) / 2, (q[0][2]+q[-1][2]) / 2)return not (point_in_tri(p, t1[0], t1[1], t1[2]) or point_in_tri(p, t2[0], t2[1], t2[2]))def gen_tri(p, t):while True:s, ok = sample([i for i in range(len(p))], 3), Truetc, ts = (p[s[0]], p[s[1]], p[s[2]]), sorted(s)for ti in t:tt = sorted(ti)if (ts[0] == tt[0] and ts[1] == tt[1] and ts[2] == tt[2]) or not valid((p[ti[0]], p[ti[1]], p[ti[2]]), tc):ok = Falsebreakif ok:return sdef gen():n, p, t = randint(3, V), [], []for _ in range(n):p.append(gen_point(p))m = randint(1, min(n*(n-1)*(n-2)//6, T))for _ in range(m):t.append(gen_tri(p, t))return n, m, p, tif __name__ == '__main__':with open('in.txt', 'w') as f:for _ in range(C):n, m, v, t = gen()f.write(f'{n} {m}\n')for p in v:f.write(f'{p[0]} {p[1]} {p[2]}\n')for p in t:f.write(f'{p[0]+1} {p[1]+1} {p[2]+1}\n')f.write('0 0\n')

测试数据可视化

  首先需要将所有空间三角形投影到$xy$平面来看区域分割情况，推荐使用在线工具geogebra来画。齐次还需要观察各个分割区域实际是哪个三角面片上的那一部分才可见，这可以借助matplotlib.pyplot写python脚本展示俯视视野下的各个空间三角形。

  给一个可视化python脚本例子：

# -*- coding: utf-8 -*-import matplotlib.pyplot as plt
import numpy as npif __name__ == '__main__':vert = np.array([[69, 57, 80],[44, 92, 31],[82, 23, 39],[54, 54, 97],[70, 100, 71],[59, 73, 92],[34, 47, 88],[10, 93, 74]])faces = [[vert[3], vert[0], vert[2]], [vert[6], vert[0], vert[2]], [vert[3], vert[6], vert[0]]]colors = ['r', 'g', 'b']ax = plt.figure().add_subplot(111, projection='3d')ax.set_xlabel('X')ax.set_ylabel('Y')ax.set_zlabel('Z')ax.view_init(elev = 90, azim = -90)for i, face in enumerate(faces):x, y, z = zip(*face)ax.plot_trisurf(np.array(x), np.array(y), np.array(z))plt.show()

测试数据

  测试数据生成和可视化脚本以及测试数据打包下载。

AC 代码

#include <iostream>
#include <iomanip>
#include <algorithm>
#include <cmath>
using namespace std;#define eps 1e-9
#define V 302
#define T 1005
int x[V], y[V], z[V], idx[T][3], a[T], b[T], c[T], d[T], e[T], q[T], v, t, m, n, kase = 0; double f[T], g[3*T*T];
struct seg {double y1, y2;bool operator< (const seg& rhs) const {return y1 + y2 < rhs.y1 + rhs.y2;}
} s[T<<1];void calc_tri_params(int (&i)[3], int j) {cin >> i[0] >> i[1] >> i[2]; --i[0]; --i[1]; --i[2];int x0 = x[i[0]], y0 = y[i[0]], z0 = z[i[0]], x1 = x[i[1]] - x0, y1 = y[i[1]] - y0, z1 = z[i[1]] - z0,x2 = x[i[2]] - x0, y2 = y[i[2]] - y0, z2 = z[i[2]] - z0;a[j] = y1*z2 - y2*z1; b[j] = z1*x2-z2*x1; c[j] = x1*y2-x2*y1;if (c[j] == 0) return;d[j] = a[j]*x0+b[j]*y0+c[j]*z0; e[m++] = j; f[j] = sqrt(a[j]*a[j] + b[j]*b[j] + c[j]*c[j]) / abs(c[j]);
}void seg_intersect(int x1, int y1, int x2, int y2, int x3, int y3, int x4, int y4) {int c = (x4-x3)*(y2-y1) - (y4-y3)*(x2-x1), t1 = (y4-y3)*(x1-x3)-(x4-x3)*(y1-y3), t2 = (y2-y1)*(x1-x3)-(x2-x1)*(y1-y3);if (c == 0 || (c > 0 && (t1 <= 0 || t1 >= c || t2 <= 0 || t2 >= c))|| (c < 0 && (t1 >= 0 || t1 <= c || t2 >= 0 || t2 <= c))) return;g[n++] = x1 + (x2-x1)*t1 / double(c);
}void tri_intersect(int i, int j) {int t1[][2] = {{x[idx[i][0]], y[idx[i][0]]}, {x[idx[i][1]], y[idx[i][1]]}, {x[idx[i][2]], y[idx[i][2]]}},t2[][2] = {{x[idx[j][0]], y[idx[j][0]]}, {x[idx[j][1]], y[idx[j][1]]}, {x[idx[j][2]], y[idx[j][2]]}};for (int k=0; k<3; ++k) for (int n = k<2 ? k+1 : 0, m=0; m<3; ++m)seg_intersect(t1[k][0], t1[k][1], t1[n][0], t1[n][1], t2[m][0], t2[m][1], t2[m<2 ? m+1:0][0], t2[m<2 ? m+1:0][1]);
}bool cmp(int i, int j) {int x1 = max(x[idx[i][0]], max(x[idx[i][1]], x[idx[i][2]])), x2 = max(x[idx[j][0]], max(x[idx[j][1]], x[idx[j][2]]));return x1 < x2;
}void intersect(int i, double x0, double x1) {int t[][2] = {{x[idx[i][0]], y[idx[i][0]]}, {x[idx[i][1]], y[idx[i][1]]}, {x[idx[i][2]], y[idx[i][2]]}};double y0[2], y1[2]; int c0 = 0, c1 = 0;for (int j=0; j<3; ++j) {int k = j<2 ? j+1 : 0;if (t[j][0] == t[k][0]) {if (abs(t[j][0]-x0) <= eps) y0[0] = t[j][1], y0[1] = t[k][1], c0 = 2;if (abs(t[j][0]-x1) <= eps) y1[0] = t[j][1], y1[1] = t[k][1], c1 = 2;} else {int a = min(t[j][0], t[k][0]), b = max(t[j][0], t[k][0]);if (c0 == 2 && abs(y0[0]-y0[1]) < eps) c0 = 1;if (c1 == 2 && abs(y1[0]-y1[1]) < eps) c1 = 1;if (c0 < 2 && x0 >= a-eps && x0 < b+eps)y0[c0++] = ((t[k][0] - x0) * t[j][1] + (x0 - t[j][0]) * t[k][1]) / (t[k][0] - t[j][0]);if (c1 < 2 && x1 >= a-eps && x1 < b+eps)y1[c1++] = ((t[k][0] - x1) * t[j][1] + (x1 - t[j][0]) * t[k][1]) / (t[k][0] - t[j][0]);}}s[v++] = {min(y0[0], y0[1]), min(y1[0], y1[1])}; s[v++] = {max(y0[0], y0[1]), max(y1[0], y1[1])};
}double a2(double x1, double y1, double x2, double y2, double x3, double y3) {return abs((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1));
}bool point_int_tri(double xp, double yp, int i) {int x0 = x[idx[i][0]], y0 = y[idx[i][0]], x1 = x[idx[i][1]], y1 = y[idx[i][1]], x2 = x[idx[i][2]], y2 = y[idx[i][2]];double s = a2(x0, y0, x1, y1, x2, y2);return abs(a2(x0, y0, x1, y1, xp, yp) + a2(x0, y0, x2, y2, xp, yp) + a2(x1, y1, x2, y2, xp, yp) - s) < eps*s;
}void solve() {for (int i=n=0; i<v; ++i) cin >> x[i] >> y[i] >> z[i], g[n++] = x[i];for (int i=m=0; i<t; ++i) calc_tri_params(idx[i], i);for (int i=0; i<m; ++i) for (int j=i+1; j<m; ++j) tri_intersect(e[i], e[j]);sort(e, e+m, cmp); sort(g, g+n); n = unique(g, g+n) - g;double ans = 0.;for (int i=1, j=0; i<n; ++i) {double x0 = g[i-1], x1 = g[i]; v = t = 0;if (x1-x0 < eps) continue;while (j < m && max(x[idx[e[j]][0]], max(x[idx[e[j]][1]], x[idx[e[j]][2]])) < x1-eps) ++j;for (int k=j; k<m; ++k) {if (min(x[idx[e[k]][0]], min(x[idx[e[k]][1]], x[idx[e[k]][2]])) > x0+eps) continue;q[t++] = e[k]; intersect(e[k], x0, x1);}sort(s, s+v);for (int k=1; k<v; ++k) {double y1 = s[k-1].y1, y2 = s[k-1].y2, y3 = s[k].y1, y4 = s[k].y2, xc, yc, zc; int c;if (abs(y3+y4-y1-y2) < eps) continue;abs(y3 - y1) < eps ? (xc = x0, yc = y1, c = 1) : (xc = x0+x0, yc = y1+y3, c = 2);abs(y4 - y2) < eps ? (xc += x1, yc += y2, c += 1) : (xc += x1+x1, yc += y2+y4, c += 2);xc = xc / c; yc = yc / c; c = -1; zc = -1e20;for (int p=0; p<t; ++p) if (point_int_tri(xc, yc, q[p])) {double z = (d[q[p]] - a[q[p]]*xc - b[q[p]]*yc) / ::c[q[p]];if (z > zc) c = q[p], zc = z;}if (c >= 0) ans += (x1-x0)*(y3+y4-y1-y2)*f[c];}}cout << "Case " << ++kase << ": " << ans / 2. << endl << endl;
}int main() {cout << fixed << setprecision(2);while (cin >> v >> t && v) solve();return 0;
}