求三维凸包,其中求体积的部分值得借鉴
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
#define MAXN 10005
const double eps = 1e-9;
int n, fc_cnt, leg[MAXN][MAXN];
double mx_ff, dis_ff;
inline int sig(double x)
{
return (x>eps) - (x<-eps);
}
struct node
{
double x, y, z;
node(double a = 0.0, double b = 0.0, double c = 0.0):x(a),y(b),z(c) {}
node operator - (node& a) const
{
return node(x-a.x, y-a.y, z-a.z);
}
node operator + (node& a) const
{
return node(x+a.x, y+a.y, z+a.z);
}
node operator ^ (node a) const // this * a * sin($)
{
return node(y*a.z-a.y*z, z*a.x-a.z*x, x*a.y-a.x*y);
}
double operator * (node a) const // this * a * cin($)
{
return x*a.x + y*a.y + z*a.z;
}
node cross(node a, node b)
{
return a - *this^b - *this;
}
double L()
{
return sqrt(x*x + y*y + z*z);
}
double vv6(node a, node b)
{
return (a^b)*(*this);
}
} P[MAXN];
struct FACE
{
int a, b, c;
int isOK;
FACE(int x = 0, int y = 0, int z = 0, int ik = 0):a(x),b(y),c(z), isOK(ik){}
int visit(node mm)
{
return sig(P[a].cross(P[b], P[c])*(mm - P[a])) > 0;
}
}facc[MAXN*8];
void deal(int o, int a, int b);
void dfs(int f, int o)
{
facc[f].isOK = 0;
dis_ff -= P[facc[f].a].vv6(P[facc[f].b], P[facc[f].c]);
deal(o, facc[f].b, facc[f].a);
deal(o, facc[f].a, facc[f].c);
deal(o, facc[f].c, facc[f].b);
}
void deal(int o, int a, int b)
{
int i = leg[a][b];
if (facc[i].isOK)
{
if (facc[i].visit(P[o]))
dfs(i, o);
else
{
dis_ff += P[b].vv6(P[a], P[o]);
FACE add(b, a, o, 1);
leg[b][a] = leg[a][o] = leg[o][b] = fc_cnt;
facc[fc_cnt++] = add;
}
}
}
void _update(int idx)
{
for (int i = 0; i< fc_cnt; i++)
{
if (facc[i].isOK && facc[i].visit(P[idx]))
{
dfs(i, idx);
break;
}
}
}
double getV()
{
double v = 0;
for (int i = 0; i<fc_cnt; i++)
{
if (facc[i].isOK)
{
v += P[facc[i].a].vv6(P[facc[i].b], P[facc[i].c]);
}
}
return v/6;
}
void solve()
{
mx_ff = dis_ff = 0;
int kg = 1, flg = 1;
int i, j;
for (i = 1; i< n; i++)
{
if (flg == 1)
{
flg += sig((P[0]-P[i]).L());
if (flg > 1)
swap(P[1], P[i]);
}
else if (flg == 2)
{
flg += sig((P[0].cross(P[1], P[i])).L());
if (flg > 2)
swap(P[i], P[2]);
}
else if (flg == 3)
{
flg += sig(P[0].cross(P[1], P[2]) * (P[i]-P[0])) != 0;
if (flg > 3)
{
swap(P[3], P[i]);
fc_cnt = 0;
for (j = 0; j< 4; j++)
{
FACE add((j+1)%4, (j+2)%4, (j+3)%4, 1);
if (add.visit(P[j]))
swap(add.b, add.c);
leg[add.a][add.b] = leg[add.b][add.c] = leg[add.c][add.a] = fc_cnt;
facc[fc_cnt++] = add;
}
for (j = 4; j<= i; j++)
{
_update(j);
}
kg = i+1;
mx_ff = getV();
}
}
else
{
dis_ff = 0;
_update(i);
dis_ff /= 6;
if (sig(dis_ff - mx_ff) > 0)
{
mx_ff = dis_ff;
kg = i+1;
}
}
}
printf("%d %.2lf\n", kg, mx_ff);
}
int main()
{
#ifndef ONLINE_JUDGE
freopen("in.txt", "r", stdin);
#endif
int cs = 0;
while (scanf("%d", &n) != EOF)
{
for (int i = 0; i< n; ++i)
{
scanf("%lf%lf%lf", &P[i].x, &P[i].y, &P[i].z);
}
printf("Case #%d:\n", ++cs);
solve();
}
return 0;
}