文档章节

空间直线与球面相交算法

o
 osc_nfjwhlc1
发布于 07/06 11:06
字数 1214
阅读 36
收藏 0
c++

行业解决方案、产品招募中!想赚钱就来传!>>>

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++实现如下:

#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("%lf\t%lf\t%lf\n", it.x, it.y, it.z);
	}	
}

最终运行的结果:
imglink1

再次注意,我这里是把线段当成直线判断的,如果希望判断整个直线与球面的交点,应该略去最后的关于\(t\)是否在0到1之间的判断,此时应该会有两个交点。

3. 参考

  1. 空间直线同球体交点求解
o
粉丝 0
博文 55
码字总数 0
作品 0
私信 提问
加载中
请先登录后再评论。
【opencv】图形的绘制

1.矩形图像的绘制: 原函数:void cvRectangle(CvArr* img, CvPoint pt1, CvPoint pt2, CvScalar color, int thickness=1, int line_type=8,int shift=0) img就是需要绘制的图像 pt1 and pt......

其实我是兔子
2014/10/08
1.1K
1
Javascript图元绘制库--ternlight

基于HTML CANVAS API的Javascript库,提供在HTML页面上绘制图元——如流程图的能力。 目前已支持简单的矩形图元和图元间的连线(直线、直角连线两种),拖拽图元等能力。 该javascript librar...

fancimage1
2013/02/07
6.2K
1
硬实时操作系统--Raw OS

Raw-OS 起飞于2012年,Raw-OS志在制作中国人自己的最优秀硬实时操作系统。 Raw-OS 操作系统特性 内核最大关中断时间无限接近0us, s3c2440系统最大关中断时间实测0.8us。 支持idle任务级别的事...

jorya_txj
2013/03/19
6.1K
1
游戏引擎--DarkGDK

Dark游戏开发工具包是一个完整的游戏引擎技术利用最新DirectX 9.0。 微软公司制作的编游戏的链接库工具,专门配合Visual C++ 2008 Express 和 DirextX 9.0 SDK,可以编辑制作3D,2D游戏,制作...

匿名
2013/04/01
2.2K
0
磁盘空间统计工具--Disk Inventory X

Disk Inventory X 是运行于 Mac OS X 10.3 (及以上)的磁盘空间统计工具。它通过 "treemaps" 的特殊方式显示了文件及文件夹的占用空间情况。 如果你经常想知道你的磁盘空间都被什么文件占用的...

匿名
2012/11/13
1.2K
0

没有更多内容

加载失败,请刷新页面

加载更多

Python 计算 0.1+0.2≠0.3? 6 张图搞清楚原理!

点击上方 Z先生点记,加为星标 第一时间收到 Python 技术干货! “ 作者:武沛齐 出处: http://www.cnblogs.com/wupeiqi/ 本文版权归作者和博客园共有 为啥会有上述图片的现象呢?其实是由于...

zeroing1
07/28
0
0
转向边缘计算? 考虑一下

云栖号资讯:【点击查看更多行业资讯】 在这里您可以找到不同行业的第一手的上云资讯,还在等什么,快来! 数据为王,特别是在当前数据驱动业务的时代,数据思维和分析能力是决定未来成功的重...

osc_lmp76vjx
2分钟前
0
0
工信部为“新基建”安全加把“锁”

云栖号资讯:【点击查看更多行业资讯】 在这里您可以找到不同行业的第一手的上云资讯,还在等什么,快来! 8月3日,工业和信息化部发布通知,要求开展2020年网络安全技术应用试点示范工作。按...

osc_oz0d1seh
3分钟前
0
0
今天吃粽子了吗?🤔祝大家端午安康~

本文分享自微信公众号 - 电子狂人(DZKR666)。 如有侵权,请联系 support@oschina.cn 删除。 本文参与“OSC源创计划”,欢迎正在阅读的你也加入,一起分享。...

狂人V
06/25
0
0
如何计算C#中某人的年龄? - How do I calculate someone's age in C#?

问题: 给定一个代表一个人生日的DateTime ,我如何计算他们的年龄(以岁为单位)? 解决方案: 参考一: https://stackoom.com/question/9/如何计算C-中某人的年龄 参考二: https://oldbug...

技术盛宴
5分钟前
13
0

没有更多内容

加载失败,请刷新页面

加载更多

返回顶部
顶部