(转)从零实现3D图像引擎：(8)参数化直线与3D平面函数库

1. 数学分析

1) 参数化直线

还记得我们在学习向量时介绍的位移向量吗？参数化直线的原理和他相同。其实所谓参数化直线，就是通过一个参数t来表示直线上的一个线段，当t为0时，则取线段的一端的端点，当t取1时，则取到另一端的端点。如图：

核心原理就是向量加法的几何意义，Vd是加数，通过t(从0到1)去乘以这个Vd，则可以描述向量p，而p在某个t时的坐标值，就是从点P1到点P2的线段在t时的坐标值。注意这里我们使用的是图中的p = p1 + t * vd，这样t的取值是0到1，如果是p = p1 + t * v'，那t的取值就是0到|Vd|，还需要去先算|Vd|，没有必要这么做，而且0和1是特殊的数字，运算方便。

所以参数化方程为：P(x,y,z) = p0 + v*t

我们为什么要用参数化的形式表示线段呢，因为在求两个线段的交点时，只要求两个参数方程组成的方程组，并求出t1和t2，只要他们都在0到1之间，那么就可以知道他们是有交点的。比如：

线段1：

x = 1 + 2*t1;

y = 3 + 4*t1;

线段2：

x = 5 + 6*t2;

y = 7 + 8*t2;

四个方程，四个未知数，可求得t1,t2，如果他们都在0和1之间，那他们相交。

具体怎么求？忘记上一章所讲的求逆矩阵的吗？要先算出行列式，然后。。。

注意，这里我们有很多的情况还没有考虑，比如：

1. 两条直线平行，没有交点

2. 两条线段在同一直线上，但互不重叠

3. 两条线段在同一直线上，部分重叠

4. 两条线段在同一直线上，只有一点重叠

在我的实现函数中，我把这4种特殊情况都看做是不相交或者说是无意义，当他们平行或共线时都会返回回传码0。因为这时计算出来的t1和t2以及交点都是没有任何意义的。

2) 3D平面

我们如何定义一个3D平面呢，看下图：

如果你认真看图了，那么无需我多说了，只要有一个点P0，和平面的法线向量n，就可以确定一个平面了。

定义：对于空间中的任意一点P，如果线段P->P0与法线向量n垂直，则点P在平面上。

还记得向量点积吗？其最最重要的性质就是其几何意义：可以确定两个向量之间的夹角。当两个向量的点积为0时，这两个向量互相垂直。

所以：n.(P->P0) = 0

展开：<a, b, c> . <x-x0, y-y0, z-z0> = 0

最终，可以拿到平面的点-法线表示的形式：a*(x-x0) + b*(y-y0) + c(z-z0) = 0

我们可以另d = (-a*x0 - b*y0 - c*z0)

那么上面的形式可以化为：

a*x + b*y + c*z + d = 0

这个就是3D平面的通用形式，以通用形式计算两个平面之间的交线非常方便。

再扩展一下思维，就可以知道，对于任意一点<x,y,z>，a*(x-x0) + b*(y-y0) + c(z-z0)的值与0的关系，可以判断出，这个点在平面的哪一边。而这个性质则非常非常的重要，我打算先举个实例再给出答案。

还是参照上图，这是个特殊的平面，是X-Z平面，所以法线向量n为<0,1,0>，平面上的一点，我们选原点<0, 0, 0>，那么现在我们来取任意一点，比如(5, 5, 5)，那么把这些数据代入a*(x-x0) + b*(y-y0) + c(z-z0)，得：

1*(5-0) = 5 > 0，所以点(5, 5, 5)在平面的正方向。

再取点(5, -5, 5)，则1*(-5-0) = -5 < 0 ，所以点(5, -5, 5)在该平面的负方向。

因为这个概念经常要用到，所以在举了一点的例子之后我想再通过理论证明一下，请看下图：

这幅图是上图的平面版，其实任意一个3D点与平面的关系都可以找到这样一个角度来看。这里的X轴就是X-Z平面的水平视角中的线，Y轴的正半轴就是这个平面的法线向量，这个图左上提示了大家，u.v = |u| * |v| * Cos(theta)。左下角则根据向量范数的求法，给出了u.v是否大于0，取决于Cos(theta)的结果。图上4种颜色的点，标出了任意一点可能在平面附近的位置。还记得余弦函数线吗？不多说了。到此我们已经证明了在任意3D卦限的某点，与某3D平面，计算u.v > 0 则在平面正方向，u.v < 0 则在负方向。

3) 参数化直线与3D平面的交点

其实就是解方程组，把直线方程和3D平面方程放在一起进行简化。

参数化直线：

x = x0 + vx * t

y = y0 + vy * t

z = z0 + vz * t

其中x0, y0, z0为参数化直线上一直的一点

3D平面：

a*x + b*y + c*z + d = 0

a = nx

b = ny

c = nz

d = -nx*x0' -ny*y0' - nz*z0'

其中x0', y0', z0'是3D平面上已知的一点，所以加了'来区分直线上已知的一点

将这几个方程代入做整理变换可以得到t的表示形式：

t = -(a*x0 + b*y0 + c*z0 + d) / (a*vx + b*vy + c*vz)

通过这个t就可以知道线段与平面是否存在交点了，不用我再说了吧，看t是不是在0和1之间即可。

这几个方程也可以表示出x, y, z，即交点的坐标：

x = x0 + vx*t

y = y0 + vy*t

z = z0 + vz*t

哈哈，其实不就是把t代回原线段的形式了么。。。

2. 代码实现

1) 结构定义

首先我们定义表示参数化直线和3D平面的结构：

typedef struct PARAMLINE2D_TYPE // 参数化2D直线
{
	POINT2D p0;
	POINT2D p1;
	VECTOR2D v;
} PARAMLINE2D, *PARAMLINE2D_PTR;

typedef struct PARAMLINE3D_TYPE // 参数化3D直线
{
	POINT3D p0;
	POINT3D p1;
	VECTOR3D v;
} PARAMLINE3D, *PARAMLINE3D_PTR;

typedef struct PLANE3D_TYPE // 3D平面
{
	POINT3D p0;
	VECTOR3D n;
} PLANE3D, *PLANE3D_PTR;

2) 常用函数实现

void _CPPYIN_Math::ParamLineCreate(POINT2D_PTR p0, POINT2D_PTR p1, PARAMLINE2D_PTR p)
{
	p->p0.x = p0->x;
	p->p0.y = p0->y;
	p->p1.x = p1->x;
	p->p1.y = p1->y;
	p->v.x = p1->x - p0->x;
	p->v.y = p1->y - p0->y;
}

void _CPPYIN_Math::ParamLineCreate(POINT3D_PTR p0, POINT3D_PTR p1, PARAMLINE3D_PTR p)
{
	p->p0.x = p0->x;
	p->p0.y = p0->y;
	p->p0.z = p0->z;
	p->p1.x = p1->x;
	p->p1.y = p1->y;
	p->p1.z = p1->z;
	p->v.x = p1->x - p0->x;
	p->v.y = p1->y - p0->y;
	p->v.z = p1->z - p0->z;
}

void _CPPYIN_Math::ParamLineGetPoint(PARAMLINE2D_PTR p, double t, POINT2D_PTR pt)
{
	pt->x = p->p0.x + p->v.x * t;
	pt->y = p->p0.y + p->v.y * t;
}

void _CPPYIN_Math::ParamLineGetPoint(PARAMLINE3D_PTR p, double t, POINT3D_PTR pt)
{
	pt->x = p->p0.x + p->v.x * t;
	pt->y = p->p0.y + p->v.y * t;
	pt->z = p->p0.z + p->v.z * t;
}

int _CPPYIN_Math::ParamLineIntersect(PARAMLINE2D_PTR p1, PARAMLINE2D_PTR p2, double *t1, double *t2) // 0 平行或共线 1 线段相交 2 线段不相交
{
	double det_p1p2 = p1->v.x * p2->v.y - p1->v.y * p2->v.x;
	if (abs(det_p1p2) <= EPSILON)
		return 0;

	*t1 = (p2->v.x*(p1->p0.y - p2->p0.y) - p2->v.y*(p1->p0.x - p2->p0.x)) /det_p1p2;
	*t2 = (p1->v.x*(p1->p0.y - p2->p0.y) - p1->v.y*(p1->p0.x - p2->p0.x)) /det_p1p2;

	if ((*t1>=0) && (*t1<=1) && (*t2>=0) && (*t2<=1))
		return 1;
	else
		return 2;
}

int _CPPYIN_Math::ParamLineIntersect(PARAMLINE2D_PTR p1, PARAMLINE2D_PTR p2, POINT2D_PTR pt) // 0 平行或共线 1 线段相交 2 线段不相交
{
	double t1, t2, det_p1p2 = (p1->v.x*p2->v.y - p1->v.y*p2->v.x);

	if (abs(det_p1p2) <= EPSILON)
		return 0;

	t1 = (p2->v.x*(p1->p0.y - p2->p0.y) - p2->v.y*(p1->p0.x - p2->p0.x))/det_p1p2;
	t2 = (p1->v.x*(p1->p0.y - p2->p0.y) - p1->v.y*(p1->p0.x - p2->p0.x))/det_p1p2;

	pt->x = p1->p0.x + p1->v.x*t1;
	pt->y = p1->p0.y + p1->v.y*t1;

	if ((t1>=0) && (t1<=1) && (t2>=0) && (t2<=1))
		return 1;
	else
		return 2;
}

void _CPPYIN_Math::PlaneCreate(PLANE3D_PTR plane, POINT3D_PTR p0, VECTOR3D_PTR normal, int normalize)
{
	plane->p0.x = p0->x;
	plane->p0.y = p0->y;
	plane->p0.z = p0->z;

	if (normalize)
	{
		VECTOR3D n;
		VectorNormalize(normal, &n);
		plane->n.x = n.x;
		plane->n.y = n.y;
		plane->n.z = n.z;
	}
	else
	{
		plane->n.x = normal->x;
		plane->n.y = normal->y;
		plane->n.z = normal->z;
	}
}

double _CPPYIN_Math::PlaneWithPoint(POINT3D_PTR pt, PLANE3D_PTR plane)
{
	double hs = plane->n.x*(pt->x - plane->p0.x) + 
		plane->n.y*(pt->y - plane->p0.y) +
		plane->n.z*(pt->z - plane->p0.z);

	return hs;
}

int _CPPYIN_Math::PlaneParamLineInterset(PARAMLINE3D_PTR pline, PLANE3D_PTR plane, double *t, POINT3D_PTR pt) // 0 不相交 1 相交 2 相交在延长线上 3 线段在平面上
{
	// 求直线与平面法线向量的点积
	double plane_dot_line = VectorDot(&pline->v, &plane->n);

	if (abs(plane_dot_line) <= EPSILON) // 点积为0, 直线与面法线垂直，要么平行于平面，要不就在平面上，下面拿直线上一点测试
	{
		if (abs(PlaneWithPoint(&pline->p0, plane)) <= EPSILON)
			return 3;
		else
			return 0;
	}

	// 求t
	*t = -(plane->n.x*pline->p0.x + 
		   plane->n.y*pline->p0.y + 
		   plane->n.z*pline->p0.z -
		   plane->n.x*plane->p0.x - 
		   plane->n.y*plane->p0.y - 
		   plane->n.z*plane->p0.z) / (plane_dot_line);
   
	pt->x = pline->p0.x + pline->v.x*(*t);
	pt->y = pline->p0.y + pline->v.y*(*t);
	pt->z = pline->p0.z + pline->v.z*(*t);

	if (*t>=0.0 && *t<=1.0)
	   return 1;
	else
	   return 2;
}

3. 代码下载

完整项目代码下载：>>点击进入下载页<<

转自:http://blog.csdn.net/cppyin/archive/2011/02/10/6177213.aspx