HDU 4449 Building Design 三维凸包+空间坐标变换

题意：给你一个三维空间的一些点，让你求一个三维凸包，并且旋转这个凸包，使得一个面贴地，并且高度最高，高度相同的情况下，映射到地面的面积最小。

思路：求一个三维凸包，然后枚举每一个面作为底面，然后套一个模板就可以求出旋转之后的各个点坐标。。然后在求二维凸包的面积。这个题不能直接用变换之后的坐标求高度，会有精度损失，直接在原坐标下求点到面的距离即可~

WA了一天，最后发现是自己的二维凸包模板有问题！！f**k

#include <iostream>
#include <string>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <vector>
#include <algorithm>
using namespace std;
const double PI = acos(-1.0);
struct Point3 {
double x, y, z;
Point3(double x=0, double y=0, double z=0):x(x),y(y),z(z) { }
};
struct Point{
double x, y;
Point(double x = 0, double y = 0) : x(x), y(y) {}
Point operator + (const Point &p){ return Point(x + p.x, y + p.y); }
Point operator - (const Point &p){ return Point(x - p.x, y - p.y); }
};
bool operator < (const Point &a, const Point &b) { return a.x < b.x || (a.x == b.x && a.y < b.y); }
bool operator < (const Point3 &a, const Point3 &b) {
if(a.x != b.x) return a.x < b.x;
else if(a.y != b.y) return a.y < b.y;
else return a.z < b.z;
}
typedef Point3 Vector3;
const double eps = 1e-8;
const Vector3 e = Vector3(0, 0, 1);
const Point3 st = Point3(0, 0, 0);
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, double p) { return Vector3(A.x/p, A.y/p, A.z/p); }
int dcmp(double x) { if(fabs(x) < eps) return 0; else return x < 0 ? -1 : 1; }
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;}
bool operator == (const Point &a, const Point &b){
return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0;
}
double Dot(const Vector3& A, const Vector3& B) { return A.x*B.x + A.y*B.y + A.z*B.z; }
double Length(const Vector3& A) { return sqrt(Dot(A, A)); }
double Angle(const Vector3& A, const Vector3& B) { return acos(Dot(A, B) / Length(A) / Length(B)); }
Vector3 Cross(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); }
double Cross(const Point &A, const Point B) { return A.x * B.y - A.y * B.x; }
double Area2(const Point3& A, const Point3& B, const Point3& C) { return Length(Cross(B-A, C-A)); }
double Volume6(const Point3& A, const Point3& B, const Point3& C, const Point3& D) { return Dot(D-A, Cross(B-A, C-A)); }
double rand01() { return rand() / (double)RAND_MAX; }
double randeps() { return (rand01() - 0.5) * eps; }
Point3 add_noise(const Point3& p) {
return Point3(p.x + randeps(), p.y + randeps(), p.z + randeps());
}
Point3 read_point3() {
Point3 P;
scanf("%lf%lf%lf",&P.x, &P.y, &P.z);
return P;
}
struct Face {
int v[3];
Face(int a, int b, int c) { v[0] = a; v[1] = b; v[2] = c; }
Vector3 Normal(const vector<Point3>& P) const {
return Cross(P[v[1]]-P[v[0]], P[v[2]]-P[v[0]]);
}
int CanSee(const vector<Point3>& P, int i) const {
return Dot(P[i]-P[v[0]], Normal(P)) > 0;
}
};
vector<Face> CH3D(const vector<Point3>& P) {
int n = P.size();
vector<vector<int> > vis(n);
for(int i = 0; i < n; i++) vis[i].resize(n);

vector<Face> cur;
cur.push_back(Face(0, 1, 2));
cur.push_back(Face(2, 1, 0));
for(int i = 3; i < n; i++) {
vector<Face> next;
for(int j = 0; j < cur.size(); j++) {
Face& f = cur[j];
int res = f.CanSee(P, i);
if(!res) next.push_back(f);
for(int k = 0; k < 3; k++) vis[f.v[k]][f.v[(k+1)%3]] = res;
}
for(int j = 0; j < cur.size(); j++)
for(int k = 0; k < 3; k++) {
int a = cur[j].v[k], b = cur[j].v[(k+1)%3];
if(vis[a][b] != vis[b][a] && vis[a][b])
next.push_back(Face(a, b, i));
}
cur = next;
}
return cur;
}
vector<Point> ConvexHull(vector<Point>& p) {
sort(p.begin(), p.end());
p.erase(unique(p.begin(), p.end()), p.end());
int n = p.size();
int m = 0;
vector<Point> ch(n+1);
for(int i = 0; i < n; i++) {
while(m > 1 && Cross(ch[m-1]-ch[m-2], p[i]-ch[m-2]) <= 0) m--;
ch[m++] = p[i];
}
int k = m;
for(int i = n-2; i >= 0; i--) {
while(m > k && Cross(ch[m-1]-ch[m-2], p[i]-ch[m-2]) <= 0) m--;
ch[m++] = p[i];
}
if(n > 1) m--;
ch.resize(m);
return ch;
}
double ConvexPolygonArea(vector<Point> p) {
double area = 0;
int n = p.size();
for(int i = 1; i < n - 1; ++i) area += Cross(p[i] - p[0], p[i + 1] - p[0]);
area = fabs(area);
return area / 2;
}
Point3 get_point(Point3 st,Point3 ed,Point3 tp){//tp在直线st-en上的垂足
double t1 = Dot(tp - st, ed - st);
double t2 = Dot(ed - st, ed - st);
double t=t1/t2;
Point3 tt = (ed-st)*t;
Point3 ans=st + tt;
return ans;
}
Point3 rotate(Point3 st,Point3 ed,Point3 tp,double A){//将点tp绕st-ed逆时针旋转A度  从ed往st看为逆时针
Point3 root=get_point(st,ed,tp);
Point3 e=(ed-st)/Length(ed - st);
Point3 r=tp-root;
Point3 vec = Cross(e, r);
Point3 ans=r*cos(A)+vec*sin(A)+root;
return ans;
}
struct ConvexPolyhedron {
int n;
vector<Point3> P, P2, P3;
vector<Face> faces;
double hi, area;
bool read() {
if(scanf("%d", &n) != 1) return false;
if(n == 0) return false;
P.resize(n);
for(int i = 0; i < n; i++) P[i] = read_point3();
sort(P.begin(), P.end());
P.erase(unique(P.begin(), P.end()), P.end());
n = P.size();
P3.resize(n);
P2.resize(n);
if(n < 4) while(1);
for(int i = 0; i < n; i++)P2[i] = add_noise(P[i]);
hi = 0, area = 1e30;
faces = CH3D(P2);
return true;
}
double mindist(Point3 C, int i) {
Point3 p1 = P[faces[i].v[0]], p2 = P[faces[i].v[1]], p3 = P[faces[i].v[2]];
double ans = fabs(-Volume6(p1, p2, p3, C) / Area2(p1, p2, p3));
return ans;
}
void solve(){
vector<Point> point2, CH2;
int fn = faces.size();
for(int i = 0; i < fn; i++){
point2.clear(); CH2.clear();
for(int j = 0; j < n; j++)
P3[j] = P[j];
Vector3 ver = Cross(P[faces[i].v[0]] - P[faces[i].v[1]], P[faces[i].v[0]] - P[faces[i].v[2]]);//求面的法向量
Vector3 ver2 = Cross(ver, e);//求旋转轴
double ang = Angle(ver, e);
if(dcmp(ang) != 0 && dcmp(ang - PI) != 0){
for(int j = 0; j < n; j++)
P3[j] = rotate(st, ver2, P3[j], ang);
}
double cur_hi = 0;
for(int j = 0; j < n; j++)
cur_hi = max(cur_hi, mindist(P[j],i));
if(dcmp(cur_hi - hi) >= 0){
for(int j = 0; j < n; j++)
point2.push_back(Point(P3[j].x, P3[j].y));
CH2 = ConvexHull(point2);
double sum = ConvexPolygonArea(CH2);
if(dcmp(cur_hi - hi) > 0) area = sum;
else area = min(area, sum);
hi = cur_hi;
}
}
printf("%.3lf %.3lf\n", hi, area);
}
};
ConvexPolyhedron CP3;
int main()
{
while(CP3.read()) CP3.solve();
return 0;
}