目录
- 1. 原理推导
- 1.1. 直线公式
- 1.2. 求交
- 2. 具体实现
- 3. 参考
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} ]
这两个方程是无法直接在实际情况中使用的，毕竟很多时候都是直接给出经过直线的两个点。我在《已知线段上某点与起点的距离，求该点的坐标》这篇博文中论述过:
对于知道线段的起点(O)和终点(E)，显然方向向量为(D=E−O)。这时，根据射线的向量方程，线段上某一点P为
[P=O tD ]
很明显，直线的参数式方程与上篇博文中描述的其实是一个意思，起点(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 实现如下：
代码语言：javascript复制#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;
}
double 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("%lft%lft%lfn", it.x, it.y, it.z);
}
}最终运行的结果:
再次注意，我这里是把线段当成直线判断的，如果希望判断整个直线与球面的交点，应该略去最后的关于(t)是否在0到1之间的判断，此时应该会有两个交点。
3. 参考
- 空间直线同球体交点求解
