zoukankan      html  css  js  c++  java
  • [Luogu4724][模板]三维凸包(增量构造法)

    1.向量点积同二维,x1y1+x2y2+x3y3。向量叉积是行列式形式,(y1z2-z1y2,z1x2-x1z2,x1y2-y1x2)。

    2.增量构造法

      1)首先定义,一个平面由三个点唯一确定。一个平面是有方向的,它的法向量只有一个方向(即逆时针相邻两向量的叉积的方向)。

      2)初始时只有(p1,p2,p3)和(p3,p2,p1)两个平面(相当于两个方向相反的面组成了一个体积为0的凸包)

      3)每次加入一个新点时,以这个点为光源中心投影到凸包上,不能被照到的面在新凸包中仍然存在,否则不存在。

      4)将新点和明暗分界线上的边组成的面加入凸包。

    3.具体实现

      1)一个面是否能被照到的判断:点在这个面的正面,即平面上任意一点到这个点的向量与平面法向量的点积为正值。

      2)明暗分界线的判断:每次判断一个面能否被照到时,标记vis[p1][p2]表示p1p2这条边的某边是否被照到。若vis[p1][p2]!=vis[p2][p1]则说明它是分界线。

      3)三位面积计算:两个向量作三维叉积后将模长除2即可。

      

    4.精度问题

      精度是计几永恒不变的话题。

      考虑四点共面的情况,在上述做法中,加入其中三个点组成的平面后,第四个点会因为找不到这个平面而被忽略,事实上第四个点可能会增加这个面的面积。

      为了避免这种情况,我们将每个点都作一些微小的“扰动”,这样就几乎不可能有四点共面的情况。同时,这个扰动幅度不能过大,不能影响到输出要求的精度。

      但是同时,由于这种扰动,如果输入数据中有多个相同的点,在扰动幅度设的不好的情况下,会因为精度误差使所有坐标相同的点全部被认为使不同的点,从而导致凸包面数指数型增长。

      实测扰动范围在1e-8左右可以同时一定程度上有效避免这两个问题,但一定也能被卡掉。也就是说网上大部分代码是可以轻松被卡掉的。具体应该怎么做才能避免精度问题感觉十分困难。

     1 #include<cmath>
     2 #include<cstdio>
     3 #include<algorithm>
     4 #define rep(i,l,r) for (int i=(l); i<=(r); i++)
     5 using namespace std;
     6 
     7 const int N=2010;
     8 const double eps=1e-8;
     9 double ans;
    10 int n,cnt,tot;
    11 bool vis[N][N];
    12 double rs() {return (rand()/(double)RAND_MAX-0.5)*eps;}
    13 
    14 struct P{ double x,y,z; }p[N];
    15 double len(P a){ return sqrt(a.x*a.x+a.y*a.y+a.z*a.z); }
    16 void shake(P &a){ a.x+=rs(); a.y+=rs(); a.z+=rs(); }
    17 P operator -(const P &a,const P &b){ return (P){a.x-b.x,a.y-b.y,a.z-b.z}; }
    18 P operator *(const P &a,const P &b){ return (P){a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x}; }
    19 double mul(const P &a,const P &b){ return a.x*b.x+a.y*b.y+a.z*b.z; }
    20 
    21 struct F{ int v[3]; }f[N],tmp[N];
    22 P Normal(F &a){ return (p[a.v[1]]-p[a.v[0]])*(p[a.v[2]]-p[a.v[0]]); }
    23 double area(F &a){ return len(Normal(a))/2.; }
    24 bool see(F &a,P &b){ return mul(b-p[a.v[0]],Normal(a))>0; }
    25 
    26 void Convex(){
    27     f[++tot]=(F){1,2,3}; f[++tot]=(F){3,2,1};
    28     rep(i,4,n){
    29         int tot1=0;
    30         rep(j,1,tot){
    31             bool fl=see(f[j],p[i]);
    32             if (!fl) tmp[++tot1]=f[j];
    33             rep(k,0,2) vis[f[j].v[k]][f[j].v[(k+1)%3]]=fl;
    34         }
    35         rep(j,1,tot) rep(k,0,2){
    36             int x=f[j].v[k],y=f[j].v[(k+1)%3];
    37             if (vis[x][y] && !vis[y][x]) tmp[++tot1]=(F){x,y,i};
    38         }
    39         tot=tot1; rep(j,1,tot1) f[j]=tmp[j];
    40     }
    41 }
    42 
    43 int main(){
    44     scanf("%d",&n);
    45     rep(i,1,n) scanf("%lf%lf%lf",&p[i].x,&p[i].y,&p[i].z),shake(p[i]);
    46     Convex();
    47     rep(i,1,tot) ans+=area(f[i]); printf("%.3lf
    ",ans);
    48     return 0;
    49 }
  • 相关阅读:
    HDU多校第六场——HDU6638 Snowy Smile(线段树区间合并)
    java
    java
    java
    java
    java
    python
    appium
    python
    python
  • 原文地址:https://www.cnblogs.com/HocRiser/p/10268901.html
Copyright © 2011-2022 走看看