空间直线与球面相交算法
钉钉、微博极速扩容黑科技,点击观看阿里云弹性计算年度发布会!>>>
目录- 1. 原理推导 1.1. 直线公式
- 1.2. 求交
1. 原理推导
1.1. 直线公式
在严格的数学定义中,直线是无线延长,没有端点的线;射线是一端有端点,另外一段没有端点无线延长的线。但在具体的计算机几何实现中,不可能去找到这种无线延长,没有端点的线,所以这里直线的定义更加近于线段,如果线段选的够长,那么这个线段就可以认为是直线或者射线。
空间直线的数学定义是,已知直线L上一点\(M_0(x_0,y_0,c_0)\),以及直线L的方向向量\(s(m,n,p)\),那么空间直线L的方程为:
\[\frac{x-x_0}{m} = \frac{y-y_0}{n} = \frac{z-z_0}{p} \]以上是空间直线的标准式方程(点向式方程)。令上面式子的比值为\(t\),那么直线的参数式方程为:
\[\begin{cases} x = x_0 + m * t\\ y = y_0 + n * t\\ z = z_0 + p * t\\ \end{cases} \]这两个方程是无法直接在实际情况中使用的,毕竟很多时候都是直接给出经过直线的两个点。我在《已知线段上某点与起点的距离,求该点的坐标》这篇博文中论述过:
\[P=O+tD \]对于知道线段的起点\(O\)和终点\(E\),显然方向向量为\(D=E−O\)。这时,根据射线的向量方程,线段上某一点P为
很明显,直线的参数式方程与上篇博文中描述的其实是一个意思,起点\(O\)就是\(M_0(x_0,y_0,c_0)\),方向向量\(D\)就是\(s(m,n,p)\):
\[\begin{cases} x = O_x + D_x * t\\ y = O_y + D_y * t\\ z = O_z + D_z * t\\ \end{cases} \tag {1} \]并且,采取这种公式描述还有个好处,局势t的取值范围为0到1,否则就在直线的两个端点之外,也就很有可能不是你需要的点。
1.2. 求交
根据数学定义,已知球心坐标\(C(C_x, C_y, C_z)\)与球的半径R,球面的公式为:
\[(X-C_x)^2 + (Y-C_y)^2 + (Z-C_z)^2 = R^2 \tag{2} \]联立(1)(2)两式,最终会得到一个关于t的一元二次方程:
\[(O_x + D_x * t-C_x)^2 + ( O_y + D_y * t-C_y)^2 + (O_z + D_z * t-C_z)^2 = R^2 \]一元二次方程组的有无解,单个解,以及双解三种可能,这也符合空间直线与球面相交的直观认识,要么相切有一个交点,要么相交有两个交点,否则的话可能没有交点。
得到\(t\)后,将其带入到(1)式中,就得到想要的交点。不过注意t的范围一般是0到1,这是与直线给的起点位置与终点位置有关的。
推到这里就会发现原来全部都是高中数学知识,应该还做过题目来着。
2. 具体实现
具体的C++实现如下:
#include <iostream> #include <string> #include <vector> using namespace std; const double EPSILON = 0.0000000001; // 3D vector struct Vector3d { public: Vector3d() { } ~Vector3d() { } Vector3d(double dx, double dy, double dz) { x = dx; y = dy; z = dz; } // 矢量赋值 void set(double dx, double dy, double dz) { x = dx; y = dy; z = dz; } // 矢量相加 Vector3d operator + (const Vector3d& v) const { return Vector3d(x + v.x, y + v.y, z + v.z); } // 矢量相减 Vector3d operator - (const Vector3d& v) const { return Vector3d(x - v.x, y - v.y, z - v.z); } //矢量数乘 Vector3d Scalar(double c) const { return Vector3d(c*x, c*y, c*z); } // 矢量点积 double Dot(const Vector3d& v) const { return x * v.x + y * v.y + z * v.z; } // 矢量叉积 Vector3d Cross(const Vector3d& v) const { return Vector3d(y * v.z - z * v.y, z * v.x - x * v.z, x * v.y - y * v.x); } bool operator == (const Vector3d& v) const { if (abs(x - v.x) < EPSILON && abs(y - v.y) < EPSILON && abs(z - v.z) < EPSILON) { return true; } return false; } do 1000 uble x, y, z; }; //求解一元二次方程组ax*x + b*x + c = 0 void SolvingQuadratics(double a, double b, double c, vector<double>& t) { double delta = b * b - 4 * a * c; if (delta < 0) { return; } if (abs(delta) < EPSILON) { t.push_back(-b / (2 * a)); } else { t.push_back((-b + sqrt(delta)) / (2 * a)); t.push_back((-b - sqrt(delta)) / (2 * a)); } } void LineIntersectSphere(Vector3d& O, Vector3d& E, Vector3d& Center, double R, vector<Vector3d>& points) { Vector3d D = E - O; //线段方向向量 double a = (D.x * D.x) + (D.y * D.y) + (D.z * D.z); double b = (2 * D.x * (O.x - Center.x) + 2 * D.y * (O.y - Center.y) + 2 * D.z* (O.z - Center.z)); double c = ((O.x - Center.x)*(O.x - Center.x) + (O.y - Center.y) * (O.y - Center.y) + (O.z - Center.z) * (O.z - Center.z)) - R * R; vector<double> t; SolvingQuadratics(a, b, c, t); for (auto it : t) { if (it >= 0 && it <= 1) { points.push_back(O + D.Scalar(it)); } } } int main() { Vector3d O(20, 30, 40); Vector3d E(20, 20, 20); Vector3d Center(20, 20, 20); double R = 15; vector<Vector3d> points; LineIntersectSphere(O, E, Center, R, points); cout<<"该直线(线段)与球面有"<< points.size() <<"个交点"<<endl; for (auto it : points) { printf("%lf\t%lf\t%lf\n", it.x, it.y, it.z); } }
最终运行的结果:
再次注意,我这里是把线段当成直线判断的,如果希望判断整个直线与球面的交点,应该略去最后的关于\(t\)是否在0到1之间的判断,此时应该会有两个交点。
3. 参考
- 空间直线段和三角形相交算法
- 空间直线段和三角形相交算法
- ActionScript实现两直线相交弧跨越算法
- 快速检测空间三角形相交算法的代码实现(Devillers & Guigue算法)
- 线段相交算法——平面扫描(可用于空间连接查询过滤)
- 空间射线与三角面相交算法及实现
- [算法]检测空间三角形相交算法(Devillers & Guigue算法)
- 两直线是否相交算法
- 28、几何算法-线段相交、凸包、球面弧长
- 快速检测空间三角形相交算法的代码实现(Devillers & Guigue算法)
- 经典算法100道(2)-绘制余弦曲线和直线相交
- 空间2直线交点算法
- GIS算法(一)直线或线段与线段相交求交点的算法(C语言)
- 1.空间中直线交点坐标问题,2.已知球面三点求球心问题以及Matlab实现
- 判断两直线是否相交
- POJ1269 简单的计算几何判断直线相交
- 空间圆弧作图算法,基于OpenGL实现!
- 空间triangle和AABB的相交判定方法
- 【解题报告】 POJ 1556 The Doors -- 最短路问题 Dijkstra算法 + 直线相交
- 24_ElasticSearch TF&IDF算法以及向量空间模型