1. 程式人生 > >HDU 4449 Building Design(計算幾何 三維凸包 + 座標轉化 模板)

HDU 4449 Building Design(計算幾何 三維凸包 + 座標轉化 模板)

題目:

You are a notable architect. 
Recently, a company invites you to design their new building for them. It is not an easy task, because they make some strange claims to the new building: 
1.Its shape should be the same as a convex polyhedron which they have given to you, though you can rotate it arbitrarily. 
2.One face of the building should cling to the ground. (Of course! Have you ever seen a building only touch the ground by an edge, or a point, or even suspend in the air?) 
3.They want the building to be as striking as possible. We call the highest point of the building as the “peak”. The higher the peak is the more striking the building will be. 
4.If there are many designs meet all claims above, the occupied land of the building should be less. In another word, the building's vertical projection area on the ground should be as small as possible. 
Now, give you the relative positions of vertices of the building, please design the building and tell them H – the height of the peak, as well as S - the projection area of the building in your best design.

Input

There are several test cases in the input. 
Each test case begins with one integer n (1 ≤ n ≤ 50), indicating the number of vertices of the building. 
Then n lines follow. Each line contains three integers x, y, z ( -10000 ≤ x, y, z ≤ 10000), separated by spaces, indicating the relative positions of one vertex of the building. 
The input ends with n = 0.

Output

For each test case, output two float numbers H and S in one line, separated by one space. 
Please round the results to three digits after decimal point.

Sample Input

6
1 0 0
-1 0 0
0 1 0
0 -1 0
0 0 1
0 0 -1
0

Sample Output

1.155 1.732

題意:

求出所給多面體的凸包,並且旋轉這個凸包,使得一個面貼地,並且高度最高,最高高度相同的情況下,則要求對映到地面的面積最小。

AC程式碼:

#include<bits/stdc++.h>
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;
}