HDU 4273 Rescue （三维凸包，求重心到表面的最小距离）

HDU4273 测试模板！ 求三维凸包，然后找重点，然后求重点到每个表面的距离。

Rescue

Time Limit: 2000/1000 MS (Java/Others) Memory Limit: 32768/32768 K (Java/Others) Total Submission(s): 326 Accepted Submission(s): 234

Problem Description

I work at NASA outer space rescue team which needs much courage and patient. In daily life, I always receive a lot of mission, and I must complete it right now. Today, team leader announced me that there is a huge spaceship dropping anchor in the out space, and we should reach there for rescue. As a working principle, at first, we should check whether there are persons living in the spaceship. So we carry a kind of machine called life sensor which can sense the life phenomenon when the distance between the machine and the living is not farther than the sense radius. I have read the designing paper of the spaceship in advance. It has a form of a convex polyhedron, and we can assume it is isodense. For best control, control center of the whole ship is located at the center of the mass. It is sure that if someone is still alive, he will stay at the control center. It’s unfortunately that I find the door is stocked when I try to enter into the spaceship, so I can only sense the living out of the space ship. Now I have opened the machine and it’s time to set the sense radius of it. I wonder the minimal radius of the machine which can allowe me to check whether there are persons living in the spaceship.

Input

There are multiple test cases. The first line contains an integer n indicating the number of vertices of the polyhedron. (4 <= n <= 100) Each of the next n lines contains three integers xi, yi, zi, the coordinates of the polyhedron vertices (-10,000 <= xi, yi, zi <= 10,000). It guaranteed that the given points are vertices of the convex polyhedron, and the polyhedron is non-degenerate.

Output

For each test case, output a float number indicating the minimal radius of the machine. Your answer should accurate up to 0.001.

Sample Input

4 0 0 0 1 0 0 0 1 0 0 0 1 8 0 0 0 0 0 2 0 2 0 0 2 2 2 0 0 2 0 2 2 2 0 2 2 2

Sample Output

0.144 1.000

Source

2012 ACM/ICPC Asia Regional Changchun Online

/* ***************

Author        :kuangbin

Created Time  :2014/10/7 8:36:38

File Name     :E:\2014ACM\专题学习\计算几何\三维几何\HDU4273.cpp

************************************************ */

#include <stdio.h>

#include <string.h>

#include <iostream>

#include <algorithm>

#include <vector>

#include <queue>

#include <set>

#include <map>

#include <string>

#include <math.h>

#include <stdlib.h>

#include <time.h>

using namespace std;

const double eps = 1e-8;

const int MAXN = 550;

int sgn(double x){

	if(fabs(x) < eps)return 0;

	if(x < 0)return -1;

	else return 1;

}

struct Point3{

	double x,y,z;

	Point3(double _x = 0, double _y = 0, double _z = 0){

		x = _x;

		y = _y;

		z = _z;

	}

	void input(){

		scanf("%lf%lf%lf",&x,&y,&z);

	}

	bool operator ==(const Point3 &b)const{

		return sgn(x-b.x) == 0 && sgn(y-b.y) == 0 && sgn(z-b.z) == 0;

	}

	double len(){

		return sqrt(x*x+y*y+z*z);

	}

	double len2(){

		return x*x+y*y+z*z;

	}

	double distance(const Point3 &b)const{

		return sqrt((x-b.x)(x-b.x)+(y-b.y)(y-b.y)+(z-b.z)*(z-b.z));

	}

	Point3 operator -(const Point3 &b)const{

		return Point3(x-b.x,y-b.y,z-b.z);

	}

	Point3 operator +(const Point3 &b)const{

		return Point3(x+b.x,y+b.y,z+b.z);

	}

	Point3 operator *(const double &k)const{

		return Point3(x*k,y*k,z*k);

	}

	Point3 operator /(const double &k)const{

		return Point3(x/k,y/k,z/k);

	}

	//点乘

	double operator *(const Point3 &b)const{

		return x*b.x + y*b.y + z*b.z;

	}

	//叉乘

	Point3 operator ^(const Point3 &b)const{

		return Point3(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);

	}

};

struct CH3D{

	struct face{

		//表示凸包一个面上的三个点的编号

		int a,b,c;

		//表示该面是否属于最终的凸包上的面

		bool ok;

	};

	//初始顶点数

	int n;

	Point3 P[MAXN];

	//凸包表面的三角形数

	int num;

	//凸包表面的三角形

	face F[8*MAXN];

	int g[MAXN][MAXN];

	//叉乘

	Point3 cross(const Point3 &a,const Point3 &b,const Point3 &c){

		return (b-a)^(c-a);

	}

	//三角形面积*2

	double area(Point3 a,Point3 b,Point3 c){

		return ((b-a)^(c-a)).len();

	}

	//四面体有向面积*6

	double volume(Point3 a,Point3 b,Point3 c,Point3 d){

		return ((b-a)^(c-a))*(d-a);

	}

	//正：点在面同向

	double dblcmp(Point3 &p,face &f){

		Point3 p1 = P[f.b] - P[f.a];

		Point3 p2 = P[f.c] - P[f.a];

		Point3 p3 = p - P[f.a];

		return (p1^p2)*p3;

	}

	void deal(int p,int a,int b){

		int f = g[a][b];

		face add;

		if(F[f].ok){

			if(dblcmp(P[p],F[f]) > eps)

				dfs(p,f);

			else {

				add.a = b;

				add.b = a;

				add.c = p;

				add.ok = true;

				g[p][b] = g[a][p] = g[b][a] = num;

				F[num++] = add;

			}

		}

	}

	//递归搜索所有应该从凸包内删除的面

	void dfs(int p,int now){

		F[now].ok = false;

		deal(p,F[now].b,F[now].a);

		deal(p,F[now].c,F[now].b);

		deal(p,F[now].a,F[now].c);

	}

	bool same(int s,int t){

		Point3 &a = P[F[s].a];

		Point3 &b = P[F[s].b];

		Point3 &c = P[F[s].c];

		return fabs(volume(a,b,c,P[F[t].a])) < eps &&

			fabs(volume(a,b,c,P[F[t].b])) < eps &&

			fabs(volume(a,b,c,P[F[t].c])) < eps;

	}

	//构建三维凸包

	void create(){

		num = 0;

		face add;

    	//***********************************
    	//此段是为了保证前四个点不共面
    	bool flag = true;
    	for(int i = 1;i < n;i++){
    		if(!(P\[0\] == P\[i\])){
    			swap(P\[1\],P\[i\]);
    			flag = false;
    			break;
    		}
    	}
    	if(flag)return;
    	flag = true;
    	for(int i = 2;i < n;i++){
    		if( ((P\[1\]-P\[0\])^(P\[i\]-P\[0\])).len() > eps ){
    			swap(P\[2\],P\[i\]);
    			flag = false;
    			break;
    		}
    	}
    	if(flag)return;
    	flag = true;
    	for(int i = 3;i < n;i++){
    		if(fabs( ((P\[1\]-P\[0\])^(P\[2\]-P\[0\]))*(P\[i\]-P\[0\]) ) > eps){
    			swap(P\[3\],P\[i\]);
    			flag = false;
    			break;
    		}
    	}
    	if(flag)return;
    	//**********************************
    
    	for(int i = 0;i < 4;i++){
    		add.a = (i+1)%4;
    		add.b = (i+2)%4;
    		add.c = (i+3)%4;
    		add.ok = true;
    		if(dblcmp(P\[i\],add) > 0)swap(add.b,add.c);
    		g\[add.a\]\[add.b\] = g\[add.b\]\[add.c\] = g\[add.c\]\[add.a\] = num;
    		F\[num++\] = add;
    	}
    	for(int i = 4;i < n;i++)
    		for(int j = 0;j < num;j++)
    			if(F\[j\].ok && dblcmp(P\[i\],F\[j\]) > eps){
    				dfs(i,j);
    				break;
    			}
    	int tmp = num;
    	num = 0;
    	for(int i = 0;i < tmp;i++)
    		if(F\[i\].ok)
    			F\[num++\] = F\[i\];
    }
    //表面积
    //测试：HDU3528
    double area(){
    	double res = 0;
    	if(n == 3){
    		Point3 p = cross(P\[0\],P\[1\],P\[2\]);
    		return p.len()/2;
    	}
    	for(int i = 0;i < num;i++)
    		res += area(P\[F\[i\].a\],P\[F\[i\].b\],P\[F\[i\].c\]);
    	return res/2.0;
    }
    double volume(){
    	double res = 0;
    	Point3 tmp = Point3(0,0,0);
    	for(int i = 0;i < num;i++)
    		res += volume(tmp,P\[F\[i\].a\],P\[F\[i\].b\],P\[F\[i\].c\]);
    	return fabs(res/6);
    }
    //表面三角形个数
    int triangle(){
    	return num;
    }
    //表面多边形个数
    //测试：HDU3662
    int polygon(){
    	int res = 0;
    	for(int i = 0;i < num;i++){
    		bool flag = true;
    		for(int j = 0;j < i;j++)
    			if(same(i,j)){
    				flag = 0;
    				break;
    			}
    		res += flag;
    	}
    	return res;
    }
    //重心
    //测试：HDU4273
    Point3 barycenter(){
    	Point3 ans = Point3(0,0,0);
    	Point3 o = Point3(0,0,0);
    	double all = 0;
    	for(int i = 0;i < num;i++){
    		double vol = volume(o,P\[F\[i\].a\],P\[F\[i\].b\],P\[F\[i\].c\]);
    		ans = ans + (((o+P\[F\[i\].a\]+P\[F\[i\].b\]+P\[F\[i\].c\])/4.0)*vol);
    		all += vol;
    	}
    	ans = ans/all;
    	return ans;
    }
    //点到面的距离
    //测试：HDU4273
    double ptoface(Point3 p,int i){
    	double tmp1 = fabs(volume(P\[F\[i\].a\],P\[F\[i\].b\],P\[F\[i\].c\],p));
    	double tmp2 = ((P\[F\[i\].b\]-P\[F\[i\].a\])^(P\[F\[i\].c\]-P\[F\[i\].a\])).len();
    	return tmp1/tmp2;
    }

};

CH3D hull;

int main()

{

    //freopen("in.txt","r",stdin);

    //freopen("out.txt","w",stdout);

    while(scanf("%d",&hull.n) == 1){

		for(int i = 0;i < hull.n;i++)hull.P[i].input();

		hull.create();

		Point3 p = hull.barycenter();

		double ans = 1e20;

		for(int i = 0;i < hull.num;i++)

			ans = min(ans,hull.ptoface(p,i));

		printf("%.3lf\n",ans);

	}

    return 0;

}