World Final 2010 LA 4795 Paperweight 三维几何

最新推荐文章于 2017-02-06 21:08:45 发布

奇怪的猪猪

最新推荐文章于 2017-02-06 21:08:45 发布

阅读量1.3k

点赞数

分类专栏： WorldFinal ACM_计算几何

本文链接：https://blog.csdn.net/jingqi814/article/details/26029841

版权

ACM_计算几何同时被 2 个专栏收录

74 篇文章 0 订阅

订阅专栏

WorldFinal

3 篇文章 0 订阅

订阅专栏

题目地址：LA4795点这里

总算把白书的例题学完了，再开始坑习题吧

额先wa了，找不出错误，把eps从1-10 改成1e-8 就ac了囧...

首先求出质心。

然后枚举每一个平面，最好的简短代码长度的方式就是把点放在数组里用for语句吧~虽然暴力了一点。

然后，一定要细分和桌面的接触面是什么，如果有四点共面，那么我们要考虑的就是这四个点的凸包（其实只用调用PointIntri就好了）

没有四点共面时就看投影是否离边距离>=0.2

我的代码：

#include<iostream>
#include<cmath>
#include<cstdio>
#include<iomanip>
#include<ctime>
#include<cstdlib>
#include<vector>

using namespace std;
const double eps=1e-8;

// 1e-8?

const double PI=acos(-1.0);

int dcmp(double x)
{
    return (x>eps)-(x<-eps);
}

struct Point3
{
    double x,y,z;
    Point3(double x=0,double y=0,double z=0):x(x),y(y),z(z) {}
    
};

typedef Point3 Vector3;

Vector3 operator+(Point3 A,Point3 B)
{
    return Point3(A.x+B.x,A.y+B.y,A.z+B.z);
    
}

Vector3 operator-(Point3 A,Point3 B)
{
    return Point3(A.x-B.x,A.y-B.y,A.z-B.z);
    
}

Vector3 operator*(Point3 A,double p)
{
    return Point3(A.x*p,A.y*p,A.z*p);
    
}

Vector3 operator/(Point3 A,double p)
{
    return Point3(A.x/p,A.y/p,A.z/p);
    
}

bool operator<(Point3 &A,Point3 &B)
{
    return A.x<B.x||(A.x==B.x&&A.y<B.y)||(A.x==B.x&&A.y==B.y&&A.z<B.z);
    
}
bool operator==(Point3 &A,Point3 &B)
{
    return dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)==0&&dcmp(A.z-B.z)==0;
}
ostream & operator<<(ostream & out,Vector3 A)
{
    out<<A.x<<' '<<A.y<<' '<<A.z<<endl;
    
    return out;
}

Point3 read_point3()
{
    Point3 p;
    scanf("%lf%lf%lf",&p.x,&p.y,&p.z);
    return p;
    
}


double Dot(Vector3 A,Vector3 B)
{
    return  A.x*B.x+A.y*B.y+A.z*B.z;
}

double Length(Vector3 A)
{
    return sqrt(Dot(A,A));
}

double Angle(Vector3 A,Vector3 B)
{
    return acos(Dot(A,B)/Length(A)/Length(B));
}


Vector3 Cross(Vector3 A,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 Area2(Point3 A,Point3 B,Point3 C)
{
    return Length(Cross(B-A,C-A));
}


double DistanceToPlane(const Point3 &p,const Point3 & p0,const Vector3 &n)
{
    return fabs(Dot(p-p0,n))/Length(n);
    
}

//  n不用管方向

Point3  GetPlaneProjection(Point3 p,Point3 p0,Vector3 n)
{
    return p-n*Dot(p-p0,n)/Dot(n,n);
}

//   三点式
Point3  GetPlaneProjection(Point3 p,Point3 p0,Point3 p1,Point3 p2)
{
       Vector3 n=Cross(p1-p0,p2-p0);
    
       return p-n*Dot(p-p0,n)/Dot(n,n);
}
Point3 LinePlaneIntersection(Point3 p1,Point3 p2,Point3 p0,Vector3 n)
{
    double t=Dot(n,p0-p1)/Dot(n,p2-p1);
    return p1+(p2-p1)*t;
    
}

//  点在三角形内部
bool   PointInTri(Point3 p,Point3 p0,Point3 p1,Point3 p2)
{
    double area1=Area2(p, p0, p1);
    double area2=Area2(p, p1, p2);
    double area3=Area2(p, p2, p0);
    
    return dcmp(area1+area2+area3-Area2(p0,p1,p2))==0;
}

//  线段与三角形相交 不考虑共面情形
bool  TriSegIntersection(Point3 p0,Point3 p1,Point3 p2,Point3 A,Point3 B,Point3 & p)
{
    Vector3 n=Cross(p1-p0, p2-p0);
    if(dcmp(Dot(n,B-A))==0)  return false;
    else{
        
        double t=Dot(n,p0-A)/Dot(n,B-A);
        if(dcmp(t)<0||dcmp(t-1)>0)  return false;
        else
        {
            p=A+(B-A)*t;
            return PointInTri(p,p0,p1,p2);
        }
        
    }
}

double DistanceToLine(Point3 P,Point3 A,Point3 B)
{
    Vector3 v1=B-A;
    Vector3 v2=P-A;
    return  Length(Cross(v1, v2))/Length(v1);
    
}

double  DistanceToSegment(Point3 P,Point3 A,Point3 B)
{
    if(A==B)  return Length(P-A);
    Vector3 v1=B-A,v2=P-A,v3=P-B;
    if(dcmp(Dot(v1,v2))<0) return Length(v2);
    else if(dcmp(Dot(v1,v3))>0) return Length(v3);
    else  return Length(Cross(v1,v2))/Length(v1);
    
}

double volume6(Point3 A,Point3 B,Point3 C,Point3 D)
{
    return  Dot(D-A,Cross(B-A,C-A));
}

// 0-1随机数

double rand01() { return rand()/(double)RAND_MAX;}
double randeps() {
    //srand((unsigned) time(NULL));
    return  (rand01()-0.5)*eps;}

Point3 add_noise(Point3 P)
{
    return Point3(P.x+randeps(),P.y+randeps(),P.z+randeps());
    
}


// 三维凸包  增量法

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]]);
    }
    // f是否能看见P[i]
    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]) // (a,b)是分界线，左边对P[i]可见
                    next.push_back(Face(a, b, i));
            }
        cur = next;
    }
    return cur;
    
}

double DistanceToPlane( const Point3 &p,const Point3 & p0,const  Point3 & p1,const  Point3 & p2)
{
    return fabs(volume6(p,p0,p1,p2)/Area2(p0, p1, p2));
    
}

//  四面体质心  就是四个点的平均
Point3 Centroid(const Point3 & A,const Point3 & B,const Point3 & C,const Point3 & D)
{
    return (A+B+C+D)/4.0;
}


struct  ConvexPolyhedron
{
    int n;  // 总顶点数
    vector<Point3>  P,P2;  //  P是真实顶点， P1是扰乱后的顶点
    vector<Face>  faces;
    
    void init(int n)
    {
        this->n=n;
        P.resize(n);
        P2.resize(n);
        
        for(int i=0;i<n;i++)
        {
            P[i]=read_point3();
            P2[i]=add_noise(P[i]);
        }
        
        
        faces=CH3D(P2);
        
    }
    
    // 本题需要,写在类里面
    
    
    double mindist(Point3 C)
    {
        double ans=1e30;
        for(int i=0;i<faces.size();i++)
        {
            ans=min(ans,DistanceToPlane(C,P[faces[i].v[0]], P[faces[i].v[1]], P[faces[i].v[2]]));  // 都是P这个数组
            
        }
        return ans;
    }
    
    
    Point3 centroid()
    {
        Point3 C=Point3(0,0,0);
        
        Point3 tot=Point3(0,0,0);
        
        double totv=0;
        
        for(int i=0;i<faces.size();i++)
        {
            Point3  p0=P[faces[i].v[0]];
            Point3  p1=P[faces[i].v[1]];
            Point3  p2=P[faces[i].v[2]];
            
            double v=volume6(p0, p1, p2, C);
            tot=tot+Centroid(p0,p1,p2,C)*v;
            totv+=v;
            
        }
        
        return tot/totv;
        
    }
    
};



//  self
// p0,p1在面A,B,C 的同侧
bool SameSide(Point3 p0,Point3 p1,Point3 A,Point3 B,Point3 C)
{
    Vector3 n=Cross(B-A, C-A);
    Vector3 v1=p0-A;
    Vector3 v2=p1-A;
    
    return Dot(v1, n)*Dot(v2,n)>=0;
}

// 改写读点函数

bool  read_point3(Point3 & P)
{
   if(scanf("%lf%lf%lf",&P.x,&P.y,&P.z)==3)
   {
       return 1;
   }
   else return 0;
    
}

Point3 P[6];

Point3 F;

//   求质心
Point3 Centroid(Point3 &A,Point3 & B,Point3 &C,Point3 &D)
{
    return (A+B+C+D)/4.0;
}

Point3 Centroid(Point3 &A,Point3 & B,Point3 &C,Point3 &D,Point3 &E)  //  注意A,B,C 为四面体公共平面
{
    double  v1=volume6(A, B, C, D);
    double  v2=volume6(A, B, C, E);
    v1=fabs(v1);
    v2=fabs(v2);
    
    Point3 C1=Centroid(A, B, C, D);
    Point3 C2=Centroid(A, B, C, E);
  
    
    return (C1*v1+C2*v2)/(v1+v2);

}

bool  ok(Point3 &P,Point3 &A,Point3 & B,Point3 &C)  // P为重心
{
    Point3 prj=GetPlaneProjection(P, A, B,C);
    if(PointInTri(prj, A, B, C))
    {
        double d1=DistanceToLine(prj, A, B);
        double d2=DistanceToLine(prj, B, C);
        double d3=DistanceToLine(prj, C, A);
        
        return  dcmp(d1-0.2)>=0&&dcmp(d2-0.2)>=0&&dcmp(d3-0.2)>=0;
        
    }
    
    else  return 0;
}

bool  ok(Point3 &P,Point3 &A,Point3 & B,Point3 &C,Point3 &D)  // P为重心 A,B 为公共边
{
    Point3 prj=GetPlaneProjection(P, A, B,C);
    if(PointInTri(prj, A, B, C)||PointInTri(prj, A, B, D))
    {
        double d1=DistanceToLine(prj, C, A);
        double d2=DistanceToLine(prj, C, B);
        double d3=DistanceToLine(prj, D, A);
        double d4=DistanceToLine(prj, D, B);

        return  dcmp(d1-0.2)>=0&&dcmp(d2-0.2)>=0&&dcmp(d3-0.2)>=0&&dcmp(d4-0.2)>=0;
        
    }
    
    else  return 0;
}

bool  OnSamePlane(Point3 &A,Point3 & B,Point3 &C,Point3 &D)
{
//    Vector3 n1=Cross(B-A, C-A);
//    Vector3 n2=Cross(B-A, D-A);
//    
//    Vector3 n1n2=Cross(n1, n2);
//    Point3 Zero=Point3(0,0,0);
//    
//    if(n1n2==Zero) return 1;
//    else return  0;
    
    return dcmp(volume6(A, B, C, D))==0;
}

void  check4(Point3 A,Point3 B,Point3 D,Point3 E,vector<double> &ans,Point3 &Center)  // A,B 是A，B,C 的棱
{
    if(OnSamePlane(A, B, D, E))
    {
        if(PointInTri(A, B, D, E))
        {
            if(ok(Center,B,D,E))
            {
                ans.push_back(DistanceToPlane(F, B, D, E));
                
            }
        }
        
        if(PointInTri(B, A, D, E))
        {
            if(ok(Center,A,D,E))
            {
                ans.push_back(DistanceToPlane(F, A, D, E));
                
            }
        }
        
        else
        {
            if(ok(Center,A,B,D,E))
            {
                ans.push_back(DistanceToPlane(F, A, B, D));
            }
        }
    }
    
    
}



int main()
{
    Point3  P_cin;
    int  cnt=0;
    int  cas=0;
    while(read_point3(P_cin))
    {
        P[(cnt++)%6]=P_cin;
        if(cnt%6==0)
        {
            F=P[5];
            cas++;
            double d1=1e30,d2=0;
            //
            Point3 Center=Centroid(P[0], P[1], P[2], P[3],P[4]);
            vector<double> ans;
            
            for(int i=0;i<5;i++)
                for(int j=i+1;j<5;j++)
                    for(int k=j+1;k<5;k++)
                    {
                        Point3 A=P[i],B=P[j],C=P[k];
                        
                        bool  fine=1;
                        for(int a=0;a<5;a++)
                            for(int b=a+1;b<5;b++)
                            {
                                if(!SameSide(P[a], P[b],A, B, C))
                                {
                                    fine=0;
                                    break;
                                }
                            }
                        
                        // 不允许共面
                        for(int i=0;i<5;i++)
                        {
                            if(P[i]==A||P[i]==B||P[i]==C) continue;
                            if(OnSamePlane(P[i], A, B, C))
                            {
                                fine=0;
                                break;
                            }
                        }
                        
                        if(!fine) continue;
                        else
                        {
                            if(ok(Center,A,B,C))
                            {
                                ans.push_back(DistanceToPlane(F, A, B, C));
                        
                            }
                        }
                        
                  }
            

          
            check4(P[0],P[1],P[3],P[4],ans, Center);
            check4(P[1],P[2],P[3],P[4],ans, Center);
            check4(P[2],P[0],P[3],P[4],ans, Center);
            
            
//            for(int i=0;i<ans.size();i++)
//            {
//                cout<<"ans"<<ans[i]<<endl;
//            }
            for(int i=0;i<ans.size();i++)
            {
                if(ans[i]<d1)
                {
                    d1=ans[i];
                }
                if(ans[i]>d2)
                {
                    d2=ans[i];
                }
            }
           
            printf("Case %d: %.5f %.5f\n",cas,d1,d2);
        }
    
    }
}

LiuRujia的代码：

// LA4795 Paperweight
// Rujia Liu
#include<cstdio>
#include<cmath>
using namespace std;

const double eps = 1e-8;
int dcmp(double x) {
    if(fabs(x) < eps) return 0; else return x < 0 ? -1 : 1;
}

struct Point3 {
    double x, y, z;
    Point3(double x=0, double y=0, double z=0):x(x),y(y),z(z) { }
};

typedef Point3 Vector3;

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

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

bool read_point3(Point3& p) {
    if(scanf("%lf%lf%lf", &p.x, &p.y, &p.z) != 3) return false;
    return true;
}

// 点p到平面p0-n的距离。n必须为单位向量
double DistanceToPlane(const Point3& p, const Point3& p0, const Vector3& n) {
    return fabs(Dot(p-p0, n)); // 如果不取绝对值，得到的是有向距离
}

// 点p在平面p0-n上的投影。n必须为单位向量
Point3 GetPlaneProjection(const Point3& p, const Point3& p0, const Vector3& n) {
    return p-n*Dot(p-p0, n);
}

// 点P到直线AB的距离
double DistanceToLine(const Point3& P, const Point3& A, const Point3& B) {
    Vector3 v1 = B - A, v2 = P - A;
    return Length(Cross(v1, v2)) / Length(v1);
}

// p1和p2是否在线段a-b的同侧
bool SameSide(const Point3& p1, const Point3& p2, const Point3& a, const Point3& b) {
    return dcmp(Dot(Cross(b-a, p1-a), Cross(b-a, p2-a))) >= 0;
}

// 点在三角形P0, P1, P2中
bool PointInTri(const Point3& P, const Point3& P0, const Point3& P1, const Point3& P2) {
    return SameSide(P, P0, P1, P2) && SameSide(P, P1, P0, P2) && SameSide(P, P2, P0, P1);
}

// 四面体的重心
Point3 Centroid(const Point3& A, const Point3& B, const Point3& C, const Point3& D) {
    return (A + B + C + D)/4.0;
}

#include<algorithm>
using namespace std;

// 判断P是否在三角形A, B, C中，并且到三条边的距离都至少为mindist。保证P, A, B, C共面
bool InsideWithMinDistance(const Point3& P, const Point3& A, const Point3& B, const Point3& C, double mindist) {
    if(!PointInTri(P, A, B, C)) return false;
    if(DistanceToLine(P, A, B) < mindist) return false;
    if(DistanceToLine(P, B, C) < mindist) return false;
    if(DistanceToLine(P, C, A) < mindist) return false;
    return true;
}

// 判断P是否在凸四边形ABCD（顺时针或逆时针）中，并且到四条边的距离都至少为mindist。保证P, A, B, C, D共面
bool InsideWithMinDistance(const Point3& P, const Point3& A, const Point3& B, const Point3& C, const Point3& D, double mindist) {
    if(!PointInTri(P, A, B, C)) return false;
    if(!PointInTri(P, C, D, A)) return false;
    if(DistanceToLine(P, A, B) < mindist) return false;
    if(DistanceToLine(P, B, C) < mindist) return false;
    if(DistanceToLine(P, C, D) < mindist) return false;
    if(DistanceToLine(P, D, A) < mindist) return false;
    return true;
}

int main() {
    for(int kase = 1; ; kase++) {
        Point3 P[5], F;
        for(int i = 0; i < 5; i++)
            if(!read_point3(P[i])) return 0;
        read_point3(F);
        
        // 求重心坐标
        Point3 c1 = Centroid(P[0], P[1], P[2], P[3]);
        Point3 c2 = Centroid(P[0], P[1], P[2], P[4]);
        double vol1 = fabs(Volume6(P[0], P[1], P[2], P[3])) / 6.0;
        double vol2 = fabs(Volume6(P[0], P[1], P[2], P[4])) / 6.0;
        Point3 centroid = (c1 * vol1 + c2 * vol2) / (vol1 + vol2);
        
        // 枚举放置方案
        double mindist = 1e9, maxdist = -1e9;
        for(int i = 0; i < 5; i++)
            for(int j = i+1; j < 5; j++)
                for(int k = j+1; k < 5; k++) {
                    // 找出另外两个点的下标a和b
                    int vis[5] = {0};
                    vis[i] = vis[j] = vis[k] = 1;
                    int a, b;
                    for(a = 0; a < 5; a++) if(!vis[a]) { b = 10-i-j-k-a; break; }
                    
                    // 判断a和b是否在平面i-j-k的异侧
                    int d1 = dcmp(Volume6(P[i], P[j], P[k], P[a]));
                    int d2 = dcmp(Volume6(P[i], P[j], P[k], P[b]));
                    if(d1 * d2 < 0) continue; // 是，则放置方案不合法
                    
                    Vector3 n = Cross(P[j]-P[i], P[k]-P[i]); // 法向量
                    n = n / Length(n); // 单位化
                    
                    Point3 proj = GetPlaneProjection(centroid, P[i], n); // 重心在平面i-j-k上的投影
                    bool ok = InsideWithMinDistance(proj, P[i], P[j], P[k], 0.2);
                    if(!ok) {
                        if(d1 == 0) { // i-j-k-a四点共面。i和j一定为ABC三个顶点之一，k和a是D或者E
                            if(!InsideWithMinDistance(proj, P[i], P[k], P[j], P[a], 0.2)) continue;
                        } else if(d2 == 0) { // i-j-k-b四点共面。i和j一定为ABC三个顶点之一，k和b是D或者E
                            if(!InsideWithMinDistance(proj, P[i], P[k], P[j], P[b], 0.2)) continue;
                        } else
                            continue;
                    }
                    
                    // 更新答案
                    double dist = DistanceToPlane(F, P[i], n);
                    mindist = min(mindist, dist);
                    maxdist = max(maxdist, dist);
                }
        printf("Case %d: %.5lf %.5lf\n", kase, mindist, maxdist);
    }
    return 0;
}