aoj1313 Intersection of Two Prisms

题意:有一个侧棱与z轴平行的棱柱$P_1$和一个侧棱与y轴平行的棱柱$P_2$。他们都向两端无限延伸,底面分别是包含M个顶点和N个顶点的凸多边形,其中第i个顶点的坐标分别是$(X_{1i},Y_{1i})$和$(X_{2i}, Z_{2i})$。计算这两个棱柱公共部分的体积。

在$x$轴上对两棱柱切片,可以得到切到的两棱柱的宽度,并且两者相乘是矩形。在$x$轴方向上对面积求积分即可。其中求积分可用抛物线法求近似值,即$\int_a^b f(x)dx \approx \frac{b-a}{6}(f(a)+4f(\frac{a+b}{2})+f(b))$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#include <bits/stdc++.h>
using namespace std;
const int maxm=100+5;
const int maxn=100+5;
int m,n;
double X1[maxn],Y1[maxn],X2[maxn],Z2[maxn];
int getmod(int id,int mo)
{
if(id==mo) return mo;
else return id%mo;
}
double wid(double *X,double *Y,int N,double x)
{
double x1,x2,y1,y2;
double L=1e20,R=-1e20,tmp=0;
for(int i=1;i<=N;i++) {
x1=X[i],y1=Y[i],x2=X[getmod(i+1,N)],y2=Y[getmod(i+1,N)];
if((x-x1)*(x-x2)<=0) {
tmp=y1+(y2-y1)*(x-x1)/(x2-x1);
L=min(L,tmp);
R=max(R,tmp);
}
}
return max(0.0,R-L);
}
vector<double > xs;
int main()
{
double L1,L2,R1,R2;
while(scanf("%d%d",&m,&n)==2) {
if(!m && !n) break;
xs.clear();
L1=1e20,L2=1e20,R1=-1e20,R2=-1e20;
for(int i=1;i<=m;i++) scanf("%lf%lf",&X1[i],&Y1[i]),xs.push_back(X1[i]),L1=min(L1,X1[i]),R1=max(R1,X1[i]);
for(int i=1;i<=n;i++) scanf("%lf%lf",&X2[i],&Z2[i]),xs.push_back(X2[i]),L2=min(L2,X2[i]),R2=max(R2,X2[i]);
sort(xs.begin(),xs.end());
double a,b,c,f1,f2,f3,ans=0.0;
for(int i=0;i+1<(int)xs.size();i++) {
a=xs[i],b=xs[i+1],c=(a+b)/2.0;
if(c>=L1 && c<=R1 && c>=L2 && c<=R2) {
f1=wid(X1,Y1,m,a)*wid(X2,Z2,n,a);
f2=wid(X1,Y1,m,b)*wid(X2,Z2,n,b);
f3=wid(X1,Y1,m,c)*wid(X2,Z2,n,c);
ans+=(b-a)/6*(f1+4*f3+f2);
}
}
printf("%.10f\n",ans);
}
return 0;
}