定义
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 +(Vector3 A,Vector3 B){return Vector3(A.x+B.x,A.y+B.y,A.z+B.z);}
Vector3 operator -(Point3 A,Point3 B){return Vector3(A.x-B.x,A.y-B.y,A.z-B.z);}
Vector3 operator *(Vector3 A,double p){return Vector3(A.x*p,A.y*p,A.z*p);}
Vector3 operator /(Vector3 A,double p){return Vector3(A.x/p,A.y/p,A.z/p);}
直线的表示,可以用参数方程 (点和向量)来表示,射线和线段可以看成”参数由取值范围限制的“的直线
平面的表示,用点法式 \((p_0,n)\), 其中 \(p_0\) 是平面上的一个点,向量 \(n\) 是平面的法向量。每个平面把空间分成了两个部分,我们用点法式表示其中一个半空间,是法向量所背离的哪一个。\(n\) 垂直于平面上的所有直线。平面上任意点 \(p\) 满足 \(Dot(n,p-p_0)=0\) 。设 \(p\) 的坐标为 \((x,y,z)\) ,\(p_0\) 的坐标为 \((x_0,y_0,z_0)\) ,向量 \(n\) 的坐标表示为 \((A,B,C)\) 上述等式等价于 $$A(x-x_0)+B(y-y_0)+C(z-z_0)=0$$
整理得: \(Ax+By+Cz-(Ax_0+By_0+Cz_0)=0\) 如果令 \(-(Ax_0+By_0+Cz_0)\),就得到了平面的一般式 $$Ax+By+Cz+D=0$$
当 \(Ax+By+Cz+D>0\) 上述点积大于 \(0\) ,即点 \((x,y,z)\) 在半平面空间 \((p_0,n)\) 外。换句话说 \(Ax+By+Cz+D>0\) 表示的是以一个半空间
点积
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));}
点到平面的距离
把 \(p-p_0\) 投影到向量 \(n\) 上可得 :\(p\) 到平面的有向距离为 \(Dot(p-p_0,n)/Length(n)\)。
double DistanceToPlane(Point3 p,Point3 p0,Vector3 n){return fabs(Dot(p-p0,n))/Length(n);}//点 p 到平面 p0-n 的距离,如果不取绝对值,得到的是有向距离
点在平面上的投影点
设 \(p\) 在平面 \((p_0,n)\) 上的投影为 \(p'\) 则 \(p'-p\) 平行于 \(n\) ,且 \(p'-p=dn\) 其中 \(d\) 为 \(p\) 到平面的的有向距离
Point3 GetPlaneProjection(Point3 p,Point3 p0,Vector3 n){//点在平面上的投影
double d=Dot(p-p0,n)/Length(n);
return p-n*d;
}
直线和平面的交点
设平面方程为 \(\operatorname{Dot}(n,p-p_0)=0\),过点 \(p_1\) 和 \(p_2\) 的直线的参数方程为 \(p=p_1+t(p_2-p_1)\) ,则与平面方程联立解得$$t=\operatorname{Dot}(n,p_0-p_1)/\operatorname{Dot}(n,p_2-p_1)$$
如果分母为 \(0\) 则直线和平面平行,或者直线在平面上
如果平面用一般式 \(Ax+By+Cz+D=0\) ,则联立出来的表达式为
\[t=\frac{Ax_1+By_1+Cz_1+D}{A(x_1-x_2)+B(y_1-y_2)+C(z_1-z_2)} \]叉积
三维叉积的结果是一个向量
\[v_1 \times v_2=\begin{bmatrix}y_1z_2-y_2z_1\\z_1x_2-z_2x_1\\x_1y_2-x_2y_1\end{bmatrix} \]叉积同时垂直于 \(v_1\) 和 \(v_2\) ,方向遵循右手定则。当且仅当 \(v_1\) 和 \(v_2\) 平行时,叉积为 \(0\)
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);}
通过向量叉积,可以进行一些拓展的基础问题。
过不共线的三点的平面,法向量为 \(\operatorname{Cross}(p_2-p_0,p_1-p_0)\) ,任取一个点即可得到平面的点法式。
三角形的有向面积的两倍
double Area2(Point3 A,Point3 B,Point3 C) {return Length(Cross(B-A,C-A));}
判断点是否在三角形内
先判断点是否在三角形所在平面上,然后利用简单的面积关系即可
bool TriSegIntersection(Point3 P0,Point3 P1,Point3 P2,Point3 A,Point3 B,Point3& P){//▲P0P1P2是否和线段AB相交
Vector3 n=Cross(P1-P0,P2-P0);
if(dcmp(Dot(n,B-A))==0) return false; //线段 AB 和平面 P0P1P2 平行或共面
else{
double t=Dot(n,P0-A)/Dot(n,B-A);
if(dcmp(t)<0||dcmp(t-1)>0) return false; //交点不在线段AB上
P=A+(B-A)*t; //计算交点
return PointInTri(P,P0,P1,P2); //判断交点是否在三角形 P0-P1-P2 内
}
}
点到直线/线段的距离,用面积法
double DistanceToLine(Point3 P,Point A,Point3 B){
Vector3 v1=B-A,V2=P-A;
return Length(Cross(v1,v2))/Length(v1);
}
四面体的体积
已知 \(4\) 个顶点为 \(A,B,C,D\)
\[V=\frac{1}{3}S\cdot h = \frac{1}{6}(\overrightarrow{AB}\times \overrightarrow{AC})\cdot h=\frac{1}{6}(\overrightarrow{AB}\times \overrightarrow{AC})\cdot \overrightarrow{AD} \]三维凸包
求三位凸包的求法有很多,常用的有暴力法,卷包裹法,和增量法。
初始时随机选取两个点 \(p_1,p_2\) 然后找一个不和两个点共线的点 \(p_3\) 找一个和他们不共面的点 \(p_4\) 组成初始凸包。 然后依次考虑其他店 \(p_r\)
- 如果这个点在当前的凸包内,直接忽略
- 否则找到这个点能 "看到" 的所有面,删除他们,然后把阴影边界上的所有点和 \(p_r\) 连接起来,其中每条边和 \(p_r\) 一起构成一个三角形
实现方法就是遍历所有面,判断是否可见,然后遍历所有边,判断是否在阴影边界上(即该边的两侧恰好有一个面可见)。
struct Face{
int v[3]; //v表示三个点的序号
Vector3 normal(Point3 *P) const{return Cross(P[v[1]]-P[v[0]],P[v[2]]-P[v[0]]);}
int cansee(Point3 *P,int i) const{return Dot(P[i]-P[v[0]],normal(P))>0?1:0;} //i能不能看到P这个平面
};
vector<Face> CH3D(Point3 *P,int n){
vector<Face> cur; //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[k].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;
}
如果存在凸包上多点共面,简单起见,实践中常常把输入点进行微笑扰动后在调用上述过程
double rand01() {return rand()/(double)RAND_MAX;}
double randeps() {return (rand01()-0.5)*eps;} //-eps/2到eps/2的随机数
Point3 add_noise(Point3 p){return Point3(p.x+randeps(),p.y+randeps(),p.z+randeps());}
扰动之前的点需要备份。上述代码记录的是 \(Face\) 是顶点的下标,所以可以用扰动后的点求三维图包,然后用这些下标引用原来 (未扰动的点)。此外,输入点应该先判重,否则在扰动后容易出现极小的面,这样的面容易造成麻烦。