UVA 11769 All Souls Night 的三维凸包要求的表面面积

主题链接：点击打开链接

求给定的3维坐标的凸包的表面积

三维凸包裸题。。

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <queue>
#include <algorithm>
using namespace std;
#define PR 1e-8
#define N 510
struct TPoint{
	double x, y, z;
	TPoint(){}
	TPoint(double _x, double _y, double _z):x(_x), y(_y), z(_z){}
	TPoint operator-(const TPoint p){return TPoint(x-p.x, y-p.y, z-p.z);}
	TPoint operator*(const TPoint p){return TPoint(y*p.z-z*p.y, z*p.x-x*p.z, x*p.y-y*p.x);}
	double operator^(const TPoint p){return x*p.x+y*p.y+z*p.z;}
};
struct fac{
	int a, b, c;
	bool ok;
};
struct T3dhull{
	int n;
	TPoint ply[N];
	int trianglecnt;
	fac tri[N];
	int vis[N][N];
	double dist(TPoint a){return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);}
	double area(TPoint a, TPoint b, TPoint c)
	{ return dist((b-a)*(c-a));}
	double volume(TPoint a, TPoint b, TPoint c, TPoint d)
	{ return (b-a)*(c-a)^(d-a);}
	double ptoplane(TPoint &p, fac &f)
	{
		TPoint m = ply[f.b] - ply[f.a], n = ply[f.c]-ply[f.a], t = p-ply[f.a];
		return (m*n)^t;
	}
	void deal(int p, int a, int b){
		int f = vis[a][b];
		fac add;
		if(tri[f].ok)
		{
			if((ptoplane(ply[p], tri[f])) > PR)
				dfs(p, f);
			else 
			{
				add.a = b, add.b = a, add.c = p, add.ok = 1;
				vis[p][b] = vis[a][p] = vis[b][a] = trianglecnt;
				tri[trianglecnt++] = add;
			}
		}
	}
	void dfs(int p, int cnt) {
		tri[cnt].ok = 0;
		deal(p, tri[cnt].b, tri[cnt].a);
		deal(p, tri[cnt].c, tri[cnt].b);
		deal(p, tri[cnt].a, tri[cnt].c);
	}
	bool same(int s, int e) {
		TPoint a = ply[tri[s].a], b = ply[tri[s].b], c = ply[tri[s].c];
		return fabs(volume(a,b,c,ply[tri[e].a])) < PR
			&& fabs(volume(a,b,c,ply[tri[e].b])) < PR
			&& fabs(volume(a,b,c,ply[tri[e].c])) < PR;
	}
	void construct()
	{
		int i, j;
		trianglecnt = 0;
		if(n<4) return ;
		bool tmp = true;
		for(i = 1; i < n; i++)
		{
			if((dist(ply[0]-ply[i])) > PR)
			{
				swap(ply[1], ply[i]);
				tmp = false;
				break;
			}
		}
		if(tmp)return ;
		tmp = true;
		for(i = 2; i < n; i++)
		{
			if((dist((ply[0]-ply[1])*(ply[1]-ply[i]))) > PR)
			{
				swap(ply[2], ply[i]);
				tmp = false;
				break;
			}
		}
		if(tmp) return ;
		tmp = true;
		for(i = 3; i < n; i++)
		{
			if(fabs((ply[0]-ply[1])*(ply[1]-ply[2])^(ply[0]-ply[i]))>PR)
			{
				swap(ply[3], ply[i]);
				tmp =false;
				break;
			}
		}
		if(tmp)return ;
		fac add;
		for(i = 0; i < 4; i++)
		{
			add.a = (i+1)%4, add.b = (i+2)%4, add.c = (i+3)%4, add.ok = 1;
			if((ptoplane(ply[i], add))>0)
				swap(add.b, add.c);
			vis[add.a][add.b] = vis[add.b][add.c] = vis[add.c][add.a] = trianglecnt;
			tri[trianglecnt++] = add;
		}
		for(i = 4; i < n; i++)
		{
			for(j = 0; j < trianglecnt; j++)
			{
				if(tri[j].ok && (ptoplane(ply[i], tri[j])) > PR)
				{
					dfs(i, j); break;
				}
			}
		}
		int cnt = trianglecnt;
		trianglecnt = 0;
		for(i = 0; i < cnt; i++)
		{
			if(tri[i].ok)
				tri[trianglecnt++] = tri[i];
		}
	}
	double area()
	{
		double ret = 0;
		for(int i = 0; i < trianglecnt; i++)
			ret += area(ply[tri[i].a], ply[tri[i].b], ply[tri[i].c]);
		return ret/2.0;
	}
	double volume()
	{
		TPoint p(0,0,0);
		double ret = 0;
		for(int i = 0; i < trianglecnt; i++)
			ret += volume(p, ply[tri[i].a], ply[tri[i].b], ply[tri[i].c]);
		return fabs(ret/6);
	}
}hull;

int main(){
	int Cas = 1;
	while(scanf("%d",&hull.n), hull.n){
		int i ;
		for(i = 0; i < hull.n; i++)
			scanf("%lf %lf %lf",&hull.ply[i].x, &hull.ply[i].y, &hull.ply[i].z);
		hull.construct();
		printf("Case %d: %.2lf\n", Cas++, hull.area());
	}
	return 0;
}

版权声明：本文博主原创文章，博客，未经同意不得转载。