三维凸包+点到平面距离+已知3点求平面方程
感谢:
http://blog.csdn.net/he11oworld/article/details/7912511
1 //已知3点坐标，求平面ax+by+cz+d=0; 2 void get_panel(TPoint p1,TPoint p2,TPoint p3,double &a,double &b,double &c,double &d){ 3 a = ( (p2.y-p1.y)*(p3.z-p1.z)-(p2.z-p1.z)*(p3.y-p1.y) ); 4 b = ( (p2.z-p1.z)*(p3.x-p1.x)-(p2.x-p1.x)*(p3.z-p1.z) ); 5 c = ( (p2.x-p1.x)*(p3.y-p1.y)-(p2.y-p1.y)*(p3.x-p1.x) ); 6 d = ( 0-(a*p1.x+b*p1.y+c*p1.z) ); 7 } 8 //点到平面距离 9 double dis_pt2panel(TPoint pt,double a,double b,double c,double d){ 10 return f_abs(a*pt.x+b*pt.y+c*pt.z+d)/sqrt(a*a+b*b+c*c); 11 }
1 /*==================================================*\ 2 | 3D凸包 3 | CALL: 构建凸包 = construct(); 4 \*==================================================*/ 5 #define TPN 1010 6 struct TPoint{ 7 double x,y,z; 8 void get(){scanf("%lf%lf%lf",&x,&y,&z);} 9 TPoint(){} 10 TPoint(double _x,double _y,double _z):x(_x),y(_y),z(_z){} 11 TPoint operator-(const TPoint p) {return TPoint(x-p.x,y-p.y,z-p.z);} 12 TPoint operator*(const TPoint p) {return TPoint(y*p.z-z*p.y,z*p.x-x*p.z,x*p.y-y*p.x);}//叉积 13 double operator^(const TPoint p) {return x*p.x+y*p.y+z*p.z;}//点积 14 }; 15 struct fac{ 16 int a,b,c;//凸包一个面上的三个点的编号 17 bool ok;//该面是否是最终凸包中的面 18 }; 19 struct T3dhull{ 20 int n;//初始点数 21 TPoint ply[TPN];//初始点 22 int trianglecnt;//凸包上三角形数 23 fac tri[TPN];//凸包三角形 24 int vis[TPN][TPN];//点i到点j是属于哪个面 25 void add(){ply[n++].get();} 26 double dist(TPoint a){return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);}//两点长度 27 double area(TPoint a,TPoint b,TPoint c){return dist((b-a)*(c-a));}//三角形面积*2 28 double volume(TPoint a,TPoint b,TPoint c,TPoint d){return (b-a)*(c-a)^(d-a);}//四面体有向体积*6 29 double ptoplane(TPoint &p,fac &f){//正：点在面同向 30 TPoint m=ply[f.b]-ply[f.a],n=ply[f.c]-ply[f.a],t=p-ply[f.a]; 31 return (m*n)^t; 32 } 33 void deal(int p,int a,int b){ 34 int f=vis[a][b];//与当前面(cnt)共边(ab)的那个面 35 fac add; 36 if(tri[f].ok){ 37 if((ptoplane(ply[p],tri[f]))>eps) dfs(p,f);//如果p点能看到该面f，则继续深度探索f的3条边，以便更新新的凸包面 38 else//否则因为p点只看到cnt面，看不到f面，则p点和a、b点组成一个三角形。 39 { 40 add.a=b,add.b=a,add.c=p,add.ok=1; 41 vis[p][b]=vis[a][p]=vis[b][a]=trianglecnt; 42 tri[trianglecnt++]=add; 43 } 44 } 45 } 46 void dfs(int p,int cnt){//维护凸包，如果点p在凸包外更新凸包 47 tri[cnt].ok=0;//当前面需要删除，因为它在更大的凸包里面 48 49 //下面把边反过来(先b,后a)，以便在deal()中判断与当前面(cnt)共边(ab)的那个面。即判断与当头面(cnt)相邻的3个面(它们与当前面的共边是反向的，如下图中(1)的法线朝外(即逆时针)的面130和312,它们共边13，但一个方向是13,另一个方向是31) 50 51 deal(p,tri[cnt].b,tri[cnt].a); 52 deal(p,tri[cnt].c,tri[cnt].b); 53 deal(p,tri[cnt].a,tri[cnt].c); 54 } 55 bool same(int s,int e){//判断两个面是否为同一面 56 TPoint a=ply[tri[s].a],b=ply[tri[s].b],c=ply[tri[s].c]; 57 return fabs(volume(a,b,c,ply[tri[e].a]))<eps 58 &&fabs(volume(a,b,c,ply[tri[e].b]))<eps 59 &&fabs(volume(a,b,c,ply[tri[e].c]))<eps; 60 } 61 void construct(){//构建凸包 62 int i,j; 63 trianglecnt=0; 64 if(n<4) return ; 65 bool tmp=true; 66 for(i=1;i<n;i++)//前两点不共点 67 { 68 if((dist(ply[0]-ply[i]))>eps) 69 { 70 swap(ply[1],ply[i]); tmp=false; break; 71 } 72 } 73 if(tmp) return; 74 tmp=true; 75 for(i=2;i<n;i++){//前三点不共线 76 if((dist((ply[0]-ply[1])*(ply[1]-ply[i])))>eps){ 77 swap(ply[2],ply[i]); tmp=false; break; 78 } 79 } 80 if(tmp) return ; 81 tmp=true; 82 for(i=3;i<n;i++){//前四点不共面 83 if(fabs((ply[0]-ply[1])*(ply[1]-ply[2])^(ply[0]-ply[i]))>eps){ 84 swap(ply[3],ply[i]); tmp=false; break; 85 } 86 } 87 if(tmp) return ; 88 fac add; 89 for(i=0;i<4;i++){//构建初始四面体(4个点为ply[0],ply[1],ply[2],ply[3]) 90 add.a=(i+1)%4,add.b=(i+2)%4,add.c=(i+3)%4,add.ok=1; 91 if((ptoplane(ply[i],add))>0) swap(add.b,add.c);//保证逆时针，即法向量朝外，这样新点才可看到。 92 vis[add.a][add.b]=vis[add.b][add.c]=vis[add.c][add.a]=trianglecnt;//逆向的有向边保存 93 tri[trianglecnt++]=add; 94 } 95 for(i=4;i<n;i++){//构建更新凸包 96 for(j=0;j<trianglecnt;j++){//对每个点判断是否在当前3维凸包内或外(i表示当前点,j表示当前面) 97 if(tri[j].ok&&(ptoplane(ply[i],tri[j]))>eps){//对当前凸包面进行判断，看是否点能否看到这个面 98 dfs(i,j); break;//点能看到当前面，更新凸包的面(递归，可能不止更新一个面)。当前点更新完成后break跳出循环 99 100 } 101 } 102 } 103 int cnt=trianglecnt;//这些面中有一些tri[i].ok=0，它们属于开始建立但后来因为在更大凸包内故需删除的，所以下面几行代码的作用是只保存最外层的凸包 104 trianglecnt=0; 105 for(i=0;i<cnt;i++){ 106 if(tri[i].ok) 107 tri[trianglecnt++]=tri[i]; 108 } 109 } 110 double area(){//表面积 111 double ret=0; 112 for(int i=0;i<trianglecnt;i++) 113 ret+=area(ply[tri[i].a],ply[tri[i].b],ply[tri[i].c]); 114 return ret/2; 115 } 116 double volume(){//体积 117 TPoint p(0,0,0); 118 double ret=0; 119 for(int i=0;i<trianglecnt;i++) 120 ret+=volume(p,ply[tri[i].a],ply[tri[i].b],ply[tri[i].c]); 121 return fabs(ret/6); 122 } 123 int facetri() {return trianglecnt;}//表面三角形数 124 int facepolygon(){//表面多边形数 125 int ans=0,i,j,k; 126 for(i=0;i<trianglecnt;i++){ 127 for(j=0,k=1;j<i;j++){ 128 if(same(i,j)) {k=0;break;} 129 } 130 ans+=k; 131 } 132 return ans; 133 } 134 135 };
