POJ 3528 求三维凸包表面积

也是用模板直接套的题目诶

//#pragma comment(linker, "/STACK:16777216") //for c++ Compiler
#include <stdio.h>
#include <iostream>
#include <cstring>
#include <cmath>
#include <stack>
#include <queue>
#include <vector>
#include <algorithm>
#define ll long long
#define Max(a,b) (((a) > (b)) ? (a) : (b))
#define Min(a,b) (((a) < (b)) ? (a) : (b))
#define Abs(x) (((x) > 0) ? (x) : (-(x)))
using namespace std;

const int INF = 0x3f3f3f3f;
const int MAXN = 666;
const double eps = 1e-8;

struct Point{
    double x,y,z;
    Point(){}
    Point(double xx,double yy,double zz):x(xx),y(yy),z(zz){}
    //两向量之差
    Point operator -(const Point p1){
        return Point(x-p1.x,y-p1.y,z-p1.z);
    }
    //两向量之和
    Point operator +(const Point p1){
        return Point(x+p1.x,y+p1.y,z+p1.z);
    }
    //叉乘
    Point operator *(const Point p){
        return Point(y*p.z-z*p.y,z*p.x-x*p.z,x*p.y-y*p.x);
    }
    Point operator *(double d){
        return Point(x*d,y*d,z*d);
    }
    Point operator / (double d){
        return Point(x/d,y/d,z/d);
    }
    //点乘
    double  operator ^(Point p){
        return (x*p.x+y*p.y+z*p.z);
    }
};

struct CH3D{
    struct face{
        //表示凸包一个面上的三个点的编号
        int a,b,c;
        //表示该面是否属于最终凸包上的面
        bool ok;
    };
    //初始顶点数
    int n;
    //初始顶点
    Point P[MAXN];
    //凸包表面的三角形数
    int num;
    //凸包表面的三角形
    face F[8*MAXN];
    //凸包表面的三角形
    int g[MAXN][MAXN];
    //向量长度
    double vlen(Point a){
        return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
    }
    //叉乘
    Point cross(const Point &a,const Point &b,const Point &c){
        return Point((b.y-a.y)*(c.z-a.z)-(b.z-a.z)*(c.y-a.y),
                     (b.z-a.z)*(c.x-a.x)-(b.x-a.x)*(c.z-a.z),
                     (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x)
                     );
    }
    //三角形面积*2
    double area(Point a,Point b,Point c){
        return vlen((b-a)*(c-a));
    }
    //四面体有向体积*6
    double volume(Point a,Point b,Point c,Point d){
        return (b-a)*(c-a)^(d-a);
    }
    //正:点在面同向
    double dblcmp(Point &p,face &f){
        Point m=P[f.b]-P[f.a];
        Point n=P[f.c]-P[f.a];
        Point t=p-P[f.a];
        return (m*n)^t;
    }
    void deal(int p,int a,int b){
        int f=g[a][b];//搜索与该边相邻的另一个平面
        face add;
        if(F[f].ok){
            if(dblcmp(P[p],F[f])>eps)
              dfs(p,f);
            else{
                add.a=b;
                add.b=a;
                add.c=p;//这里注意顺序,要成右手系
                add.ok=true;
                g[p][b]=g[a][p]=g[b][a]=num;
                F[num++]=add;
            }
        }
    }
    //递归搜索所有应该从凸包内删除的面
    void dfs(int p,int now){
         F[now].ok=0;
         deal(p,F[now].b,F[now].a);
         deal(p,F[now].c,F[now].b);
         deal(p,F[now].a,F[now].c);
    }
    bool same(int s,int t){
        Point &a=P[F[s].a];
        Point &b=P[F[s].b];
        Point &c=P[F[s].c];
        return fabs(volume(a,b,c,P[F[t].a]))<eps &&
               fabs(volume(a,b,c,P[F[t].b]))<eps &&
               fabs(volume(a,b,c,P[F[t].c]))<eps;
    }
    //构建三维凸包
    void create(){
        int i,j,tmp;
        face add;
        num=0;
        if(n<4)return;
    //**********************************************
        //此段是为了保证前四个点不共面
        bool flag=true;
        for(i=1;i<n;i++){
            if(vlen(P[0]-P[i])>eps){
                swap(P[1],P[i]);
                flag=false;
                break;
            }
        }
        if(flag)return;
        flag=true;
        //使前三个点不共线
        for(i=2;i<n;i++){
            if(vlen((P[0]-P[1])*(P[1]-P[i]))>eps){
                swap(P[2],P[i]);
                flag=false;
                break;
            }
        }
        if(flag)return;
        flag=true;
        //使前四个点不共面
        for(int i=3;i<n;i++){
            if(fabs((P[0]-P[1])*(P[1]-P[2])^(P[0]-P[i]))>eps){
                swap(P[3],P[i]);
                flag=false;
                break;
            }
        }
        if(flag)return;
    //*****************************************
        for(i=0;i<4;i++){
            add.a=(i+1)%4;
            add.b=(i+2)%4;
            add.c=(i+3)%4;
            add.ok=true;
            if(dblcmp(P[i],add)>0)swap(add.b,add.c);
            g[add.a][add.b]=g[add.b][add.c]=g[add.c][add.a]=num;
            F[num++]=add;
        }
        for(i=4;i<n;i++){
            for(j=0;j<num;j++){
                if(F[j].ok&&dblcmp(P[i],F[j])>eps){
                    dfs(i,j);
                    break;
                }
            }
        }
        tmp=num;
        for(i=num=0;i<tmp;i++)
          if(F[i].ok)
            F[num++]=F[i];
    }
    //表面积
    double area(){
        double res=0;
        if(n==3){
            Point p=cross(P[0],P[1],P[2]);
            res=vlen(p)/2.0;
            return res;
        }
        for(int i=0;i<num;i++)
          res+=area(P[F[i].a],P[F[i].b],P[F[i].c]);
        return res/2.0;
    }
    double volume(){
        double res=0;
        Point tmp(0,0,0);
        for(int i=0;i<num;i++)
           res+=volume(tmp,P[F[i].a],P[F[i].b],P[F[i].c]);
        return fabs(res/6.0);
    }
    //表面三角形个数
    int triangle(){
        return num;
    }
    //表面多边形个数
    int polygon(){
        int i,j,res,flag;
        for(i=res=0;i<num;i++){
            flag=1;
            for(j=0;j<i;j++)
              if(same(i,j)){
                  flag=0;
                  break;
              }
            res+=flag;
        }
        return res;
    }
    //三维凸包重心
    Point barycenter(){
        Point ans(0,0,0),o(0,0,0);
        double all=0;
        for(int i=0;i<num;i++){
            double vol=volume(o,P[F[i].a],P[F[i].b],P[F[i].c]);
            ans=ans+(o+P[F[i].a]+P[F[i].b]+P[F[i].c])/4.0*vol;
            all+=vol;
        }
        ans=ans/all;
        return ans;
    }
    //点到面的距离
    double ptoface(Point p,int i){
        return fabs(volume(P[F[i].a],P[F[i].b],P[F[i].c],p)/vlen((P[F[i].b]-P[F[i].a])*(P[F[i].c]-P[F[i].a])));
    }
};
CH3D hull;
int main(){
    while(scanf("%d",&hull.n)==1){
        for(int i=0;i<hull.n;i++){
            scanf("%lf%lf%lf",&hull.P[i].x,&hull.P[i].y,&hull.P[i].z);
        }
        hull.create();
        printf("%.3f\n",hull.area());
    }
    return 0;
}
时间: 2024-10-10 22:46:16

POJ 3528 求三维凸包表面积的相关文章

三维凸包模版 求三维凸包的表面积和体积

#include <stdio.h> #include <string.h> #include <stdlib.h> #include <math.h> #include <iostream> #include <queue> #include <algorithm> using namespace std; #define PR 1e-8 #define N 510 struct TPoint{ double x, y,

UVA1438 Asteroids(增量法求三维凸包,加权所有三棱锥质量求多面体重心)

https://www.luogu.com.cn/problem/UVA1438 题解建议参考lrj白书. 代码来自牛逼网友. 我只是存个板子. 1 #include <cstdio> 2 #include <cstring> 3 #include <cmath> 4 #include <cstdlib> 5 #include <vector> 6 #include <algorithm> 7 8 using namespace st

POJ 3528 hdu 3662 三维凸包模板题

POJ 3528题:http://poj.org/problem?id=3528 HDU 3662:http://acm.hdu.edu.cn/showproblem.php?pid=3662 一个是求三维凸包面数,一个是求三维凸包表面积,都是很裸的. 贴代码: #include<stdio.h> #include<algorithm> #include<string.h> #include<math.h> #include<stdlib.h>

POJ 2225 / ZOJ 1438 / UVA 1438 Asteroids --三维凸包,求多面体重心

题意: 两个凸多面体,可以任意摆放,最多贴着,问他们重心的最短距离. 解法: 由于给出的是凸多面体,先构出两个三维凸包,再求其重心,求重心仿照求三角形重心的方式,然后再求两个多面体的重心到每个多面体的各个面的最短距离,然后最短距离相加即为答案,因为显然贴着最优. 求三角形重心见此: http://www.cnblogs.com/whatbeg/p/4234518.html 代码:(模板借鉴网上模板) #include <iostream> #include <cstdio> #in

UVA 11769 All Souls Night 三维凸包求表面积

题目链接:点击打开链接 求给定的3维坐标的凸包的表面积 三维凸包裸题.. #include <stdio.h> #include <string.h> #include <stdlib.h> #include <math.h> #include <iostream> #include <queue> #include <algorithm> using namespace std; #define PR 1e-8 #def

POJ 3528

三维凸包 /* 增量法求凸包.选取一个四面体,同时把它各面的方向向量向外,增加一个点时,若该点与凸包上的某些面的方 向向量在同一侧,则去掉那些面,并使某些边与新增点一起连成新的凸包上的面. */ #include <iostream> #include <cstring> #include <cstdio> #include <algorithm> #include <cmath> using namespace std; const int M

POJ3528 HDU3662 三维凸包模板

POJ3528 HDU3662 第一道题 给定若干点 求凸包的表面积,第二题 给定若干点就凸包的面数. 简单说一下三维凸包的求法,首先对于4个点假设不共面,确定了唯一四面体,对于一个新的点,若它不在四面体内,为了让它进入凸包, 则对于所有凸包上的边,若边的一面是该点可以看到的而另一面看不到,则该点与该边构成的面要加入凸包. 模板代码非常清晰, #include<stdio.h> #include<algorithm> #include<string.h> #includ

hdu4266(三维凸包模板题)

/*给出三维空间中的n个顶点,求解由这n个顶点构成的凸包表面的多边形个数. 增量法求解:首先任选4个点形成的一个四面体,然后每次新加一个点,分两种情况: 1> 在凸包内,则可以跳过 2> 在凸包外,找到从这个点可以"看见"的面,删除这些面, 然后对于一边没有面的线段,和新加的这个点新建一个面,至于这个点可以看见的面, 就是求出这个面的方程(可以直接求法向量). */ #include<iostream> #include<cmath> #includ

hdu4449Building Design(三维凸包+平面旋转)

链接 看了几小时也没看懂代码表示的何意..无奈下来问问考研舍友. 还是考研舍友比较靠谱,分分钟解决了我的疑问. 可能三维的东西在纸面上真的不好表示,网上没有形象的题解,只有简单"明了"的讲解. 这题说起来很简单,求下三维凸包,枚举每一个面,进行坐标旋转,使得当前面作为xoy面时的其他坐标,然后求下投影面的凸包的面积. 为什么要旋转面而不直接算点到面的距离,是因为投影的面积没有办法算. 面旋转时是怎么旋转的,首先求得当前面的法向量p1,再求得它与向量e(0,0,1)的法向量pp,所有的点