Ultraviolet A 11769 All Souls Night three-dimensional convex hull surface area, 11769 souls
Question link: Click the open link
Evaluate the surface area of a convex hull given with a 3-dimensional coordinate
..
#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 510struct 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;}