Tetrahedrons and Spheres

Tetrahedrons and Spheres

时间: 1ms        内存:128M

描述:

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).

 

Input

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.

 

Output

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

提示:

参考答案(内存最优[876]):

#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;
}

参考答案(时间最优[588]):

#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;
}

题目和答案均来自于互联网,仅供参考,如有问题请联系管理员修改或删除。

点赞

发表评论

电子邮件地址不会被公开。必填项已用 * 标注