NEERC 2008 Aerodynamics

Link To The Problem

转化成2维的凸包

// 3D截面 + 2D凸包 

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<vector>
#include<cmath>

using namespace std;

#define nMax 110
#define sf scanf
#define pf printf
#define rep(i,n) for(int (i)=0;(i)<(n);(i)++)

#define pb push_back
#define mp make_pair

double const eps = 1e-9;
double const pi = acos(-1.0);

inline int dcmp(double x) {
	if(x>=-eps && x<=eps) return 0;
	return x > 0 ? 1 : -1;
}

/***********基础*************/
struct Point3 {
  double x, y, z;
  Point3(double x=0, double y=0, double z=0):x(x),y(y),z(z) { }
  void read() { sf("%lf%lf%lf",&x,&y,&z); }
  void out() { pf("(%.4lf %.4lf %.4lf)\n",x,y,z); }
};

typedef Point3 Vector3;

Vector3 operator + (const Vector3& A, const Vector3& B) { return Vector3(A.x+B.x, A.y+B.y, A.z+B.z); }
Vector3 operator - (const Point3& A, const Point3& B) { return Vector3(A.x-B.x, A.y-B.y, A.z-B.z); }
Vector3 operator * (const Vector3& A, double p) { return Vector3(A.x*p, A.y*p, A.z*p); }
Vector3 operator * (const Vector3& A, const Vector3& B) { return Vector3(A.y*B.z - A.z*B.y, A.z*B.x - A.x*B.z, A.x*B.y - A.y*B.x); }
Vector3 operator / (const Vector3& A, double p) { return Vector3(A.x/p, A.y/p, A.z/p); }
double  operator ^ (const Vector3& A, const Vector3& B) { return A.x*B.x + A.y*B.y + A.z*B.z; }
bool operator == (const Point3& A,const Point3& B)  { return dcmp(A.x-B.x)==0 && dcmp(A.y-B.y)==0 && dcmp(A.z-B.z)==0 ; }


struct point{
	double x,y;
	point() {};
	point(double x,double y):x(x),y(y){};
	point operator - (point v) { return point(x-v.x,y-v.y); }
	point operator + (point v) { return point(x+v.x,y+v.y); }
	point operator * (double k) { return point(x*k,y*k); }
	point operator / (double k) { return point(x/k,y/k); }
	double operator * (point v) { return x*v.y - y*v.x; }
	double operator ^ (point v) { return x*v.x + y*v.y; }

	void out() { pf("(%.4lf %.4lf)\n",x,y); }
};
int cmp(point a,point b) {
	if(dcmp(a.x-b.x)==0) return dcmp(a.y-b.y) < 0;
	return dcmp(a.x-b.x) < 0;
}
struct face {
	Point3 p;
	Vector3 n;
	face() {};
	face(Point3 p,Vector3 n):p(p),n(n){};
};

#define bug puts("Here");

Point3 GetPlaneProjection(const Point3& p, face s) {
  return p-s.n*((p-s.p)^s.n);
}
point f(Point3 p1,Point3 p2,face s){
	Vector3 v = p2-p1;
	if(dcmp(s.n^(p2-p1)) == 0){
		v = GetPlaneProjection(p1,s);
		return point(v.x,v.y);
	}
//	s.n.out();
	double t = (s.n^(s.p-p1))/(s.n^(p2-p1));
	Point3 ret = p1 + v*t;
	return point(ret.x,ret.y);
}
int oppiside(Point3 p1,Point3 p2,face s){
	return dcmp((s.n^(p1-s.p)) * (s.n^(p2-s.p))) <= 0;
}



double Area(point p[],int n) {
	point q[nMax*nMax];
	//rep(i,n) p[i].out();
	//cout << n << endl;
	sort(p,p+n,cmp);
	//bug
	int m=0;
	for(int i=0;i<n;i++) {
		while(m>1 && dcmp((q[m-1]-q[m-2])*(p[i]-q[m-1])) <= 0) m--;
		q[m++]=p[i];
	}
	int k = m;
	for(int i=n-2;i>=0;i--) {
		while(m>k && dcmp((q[m-1]-q[m-2])*(p[i]-q[m-1])) <= 0) m--;
		q[m++]=p[i];
	}
	double ans = 0.0;
	for(int i=1;i<m-1;i++) {
		ans += (q[i+1]-q[0])*(q[i]-q[0]);
	}
	ans /= 2.0;
	return fabs(ans);
}

Point3 p[nMax];
int n,z1,z2;
int m;
point q[nMax*nMax];

void work(int z0) {
	face s(Point3(0,0,z0),Vector3(0,0,1));
	m = 0;
	rep(i,n) for(int j=i+1;j<n;j++) {
		if(oppiside(p[i],p[j],s)){
			q[m++]=f(p[i],p[j],s);
		}
	}
	double area = Area(q,m);
	pf("%.9lf\n",area);
	return ;
}

	
int main() {
#ifndef ONLINE_JUDGE
	freopen("in.txt","r",stdin);
#endif
	
	int cas = 0;
	while(~sf("%d%d%d",&n,&z1,&z2)) {
		if(cas) pf("\n");
		rep(i,n) p[i].read();
		for(int i=z1;i<=z2;i++) {
			work(i) ;
		}
		cas ++;
	}

	return 0;
}