基于B-Spline的刚体运动

<script src="http://widgets.amung.us/classic.js" type="text/javascript"></script> <script type="text/javascript"> </script>

基于B-Spline的刚体运动

2009-09-02

 

之前我有一篇文章介绍了关于Bezier曲线和刚体运动的一些方面，现在再来介绍一下关于B-Spline曲线和刚体运动的一些东西。但不管是基于Bezier曲线的刚体运动还是基于B-Spline曲线的刚体运动，它们的核心都是一样的，都是要用一条轨迹曲线和一条方向曲线上的点,然后构成刚体变换矩阵M。顾名思义，轨迹曲线上的点决定了刚体的位置，方向曲线上的点决定了刚体的朝向。

 

 

 

对刚体运动不了解的话，可以参考我的一篇文章 Motion Analysis HTR，这篇文章对刚体运动有详细的描述。现在主要讲一下如何把B-Spline曲线和刚体运动关联起来，以及如何应用到OpenGL中来。先来看看曲线的定义，B-Spline曲线定义为

 

 

 

其中，n为最后一个控制点的索引，如果有5个控制点，最后一个索引就是4。为p次第i项B-Spline曲线基函数，为第i个控制点。对B-Spline曲线不了解的话，可以参考我的另一篇文章，Intorduction to B-Spline Curve，这篇文章对B-Spline进行了详细的介绍。当然在实际应用中，很多时候需要Nurbs，(Non-Uniform Rational B-Spline)，Nurbs没有什么神秘的地方，它也是B-Spline，只不过是knot间隔不相等的有理B-Spline。关于有理B-Spline这里只给出公式

 

 

上式中的w为权重值。

 

下面来看看在程序中的应用，这里以OpenGL为例。要画出一条曲线，首先要建立一条曲线。我们可以编写一个曲线类，比如叫做Nurbs_Curve，代码如下。

 

// _p：曲线的次数 // _n：最后一个控制点的索引 // _P：控制点数组的首地址 // _U：knot数组的首地址 Nurbs_Curve(const int _p, const int _n, const GVector *_P, const double *_U); //这样我们就可以定义一个曲线了。 //控制点数组 GVector P[] = { GVector(4, -15.0, 0.0,0.0, 1.0), GVector(4, -4, 10.0, 0.0, 1.0), GVector(4, 4.0,10.0, 0.0, 1.0), GVector(4, 15.0,0.0,0.0,1.0) }; //knot 数组 double U[] = {0.0, 0.0, 0.0, 0.0, 1.0,1.0, 1.0, 1.0}; //生成曲线 Nurbs_Curve curve(3, 3, P, U);

 

上面代码中，GVector是一个向量内，GVector(4,…)表示向量是4维向量，后面的参数就是对应的每一个分量。定义了曲线后，要画出这条曲线，就要得到曲线上的点。由曲线的方程我们知道，在knot区间上的每一个参数u都对应 一个曲线上的点，一般地，u在0和1之间。有了曲线上的点，在OpenGL中用线连接起来就可以了。

 

glLineWidth(2.0); glBegin(GL_LINE_STRIP); glColor3f(1.0,1.0,0.0); for(int i = 0; i<100; i++) { double u = curve.min_of_domain + (curve.max_of_domain-curve.min_of_domain)*i/99; GVector v = curve.Point_On_Curve(u); glVertex3f(v[0],v[1],v[2]); } glEnd();   

 

上面代码中，使用GL_LINE_STRIP是为了让每个点都可以用线连接起来，循环中100表示计算曲线上100个点。

 double u = curve.min_of_domain + (curve.max_of_domain-curve.min_of_domain)*i/99;这句代码保证参数在曲线的knot区间之内。 max_of_domain和 min_of_domain分别表示曲线knot区间的最大值和最小值。而 Point_On_Curve()是曲线内的一个成员函数，参数u是个double型的变量，它返回一个向量值，表示曲线上的点，现在就可以把这些曲线上的点连接起来。

 

 

下面来看看如何应用刚体运动。首先刚体运动是用矩阵乘法来实现的，所以我们要有矩阵类提

供基本的矩阵操作。如果矩阵类为GMatrix的话，在程序中就要申明一个全局的矩阵M作为刚

体变换矩阵，然后分别计算出位移变换矩阵T和旋转矩阵R，那么利用M=TR就可以求出这个

刚体变换矩阵了。先来看看只有位移没有旋转的情况。那么

 

 

 

这个全局变量M就实时保存了每一时刻的曲线上的位置。我们可以把这个计算变换矩阵M的过程放到一个单独的函数中完成。

 

//全局刚体变换矩阵M GMatrix M; void Motion() { //t随时间而变化,所以用静态变量,用全局变量也可以 static double t = curve.min_of_domain; //获得曲线上的点 GVector v = curve.Point_On_Curve(t); //构建变换矩阵 double e[] = { 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, v[0], v[1], v[2], 1.0 } ; GMatrix T(4, 4, e); M = T; t += 0.003; if( t> NC->max_of_domain ) t = NC->min_of_domain; }   

 

代码中,首先用一个静态变量t来作为曲线的参数u,所以要保证这个t在曲线的knot的区间之内。

然后用向量v来保存曲线上的点，把v的分量v[0],v[1],v[2]对应的x，y，z，作为矩阵的元素传入位移变换矩阵T，然后将T赋值给全局的刚体变换矩阵M。这里要注意的是，矩阵元素存放的方法和我们计算的时候不一样，这里恰好是将我们计算用的矩阵转置了一下。这是OpenGL中矩阵的结构决定的。有了变换矩阵了，如何应用它到OpenGL中呢？要在OpenGL中使用自己定义的矩阵，就不能使用OpenGL自带的glTranslatef()函数，而要使用glMultMatrixd()函数。这个函数有个参数，double型的指针。这个指针要指向我们的变换矩阵的元素数组的首地址。

//从矩阵返回矩阵元素数组首地址 double *m = form_arr(M); glPushMatrix(); //进行矩阵乘法 glMultMatrixd(m); glutSolidCube(2.0); glPopMatrix();   

代码中用m保存了从变换矩阵中返回的元素数组的地址。form_arr()这个函数可以编写为GMatrix类的友元函数，然后使用它返回GMatrix中的数组地址。有了上面的操作，我们就可以看到刚体沿着我们设计好的曲线运动。

 

 

 

 

 

上面就是在3个不同u时刻的截图，可以看到刚体的运动轨迹就是我们设计好的曲线。

 

现在来看看如何定义方向曲线和怎样应用方向曲线。如果是平面上的曲线，那么控制点是二维的。如果是空间上的曲线，那么控制点是三维的。当然控制点也可以是四维的。如果我们把这个四维向量看作为四元数的话，那么这个曲线就是方向曲线。因为单位四元数可以表示围绕任意轴的旋转。不明白什么是四元数的话，请参考其他计算机图形数的书籍。这里只给出围绕(x,y,z)轴旋转的四元数：

 

 

这样，我们在定义曲线控制点的时候，就要给出这样的四元数。所以我们要有一个四元数的类，这里为GQuater。

 

 

GQuater q1( GVector3(1.0,0.0,0.0), 0); GQuater q2( GVector3(0.0,1.0,0.0), 45); GQuater q3( GVector3(0.0,1.0,0.0), 90); GQuater q4( GVector3(0.0,0.0,1.0), 90); GVector R[] = { form_vec(q1), form_vec(q2), form_vec(q3), form_vec(q4)}; Nurbs_Curve r_curve(3,3,R,U);   

上面的代码中首先定义了四个四元数，分别表示绕x轴旋转0度，绕y轴旋转45度，绕y轴旋转90度和绕z轴旋转90度。然后把这些表示旋转的四元数用向量表示出来， form_vec()这个方法就是把四元数表示为向量，然后存放在数组R中，把R作为曲线的控制点传入曲线。这里的U和轨迹曲线的U参数是一样的。

 

现在方向曲线也有了，同样的道理可以得到方向曲线上的点来表示刚体的旋转。还有点要注意，方向曲线上的点是用向量表示的，但是我们的刚体变换矩阵M是用矩阵的乘法来实现的。所以这里要把这个向量转化为矩阵。要把方向向量表示为矩阵，首先要把它表示回四元数的形式。然后再把这个四元数表示为矩阵。四元数表示为矩阵比较复杂，这里只给出公式，要注意的是，四元数必须是单位四元数。

 

 

有了这个旋转矩阵，我们就可以通过M=TR来计算变换矩阵了。现在只要在Motion函数中添加旋转矩阵就可以了。

 

void Motion() { static double t = curve.min_of_domain; GVector v = curve.Point_On_Curve(t); //获得方向曲线上的点 GVector r = r.curve.Point_On_Curve(t); //将向量转换为四元数 GQuater q(r[3],r[0],r[1],r[2]); //将四元数转换为矩阵 GMatrix R = form_mat(q,false); double e[] = { 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, v[0], v[1], v[2], 1.0 } ; GMatrix T(4, 4, e); //计算刚体变换矩阵 M = R*T; t += 0.003; if( t> NC->max_of_domain ) t = NC->min_of_domain; }  

代码中M=R*T，这也是因为OpenGL矩阵结构的原因。注意矩阵不满足交换律的，RT和TR的结果是不一样的。

 

 

  

 

现在看看运行效果，可以明显看到不同位置的刚体方向也发生了变化。

 

这就是B-Spline和刚体运动的结合。对于曲线，不止这么简单。如果要做一个完整刚体曲线运动设计程序的话，还要涉及到曲线的一些基本算法，比如knot insertion， knot removing，knot refinement，control point reposition等等。有兴趣的朋友可以参看《The NURBS》这本书。

 

*原创文章，转载请注明出处*