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;
}