今日,做3D图形渲染,需要用到四元数,上网搜了相关的学习资料,整理如下:
(http://blog.csdn.net/pizi0475/article/details/5584617)四元数是通过使用四个数来表达方位,其目的是避免旋转过程中的万向锁问题。所以在3D中,四元数的主要用途即是用于旋转。从数学意义上讲,四元数是扩展了复数系统,它使用三个虚部i,j,k,一个四元数[x, y, z, w]定义了复数w+xi+yj+zk。我们通过OSG中的Quat类来了解四元数的基本操作及应用
1. 数据定义
使用一维数组来保存数据 double _v[4]
2. 基本操作:与向量相同,主要的基本操作都用内联函数来实现
四元数的标题乘除运算
四元数的叉乘:
四元数的叉乘即可表示为两个四元数相乘(与复数相乘方法相同),但注意叉乘只满足律,不满足交换律
且四元数乘积的模=模的乘积
设q为旋转四元数形式[cos(△/2), nsin(△/2)],n为旋转轴,单位向量;△为旋转角,执行下面的乘法可以使3D点p绕n旋转 p’=qpq-1
四元数乘法能用来连接多次旋转,与矩阵乘法效果一样。
四元数除四元数:OSG中表示为当前四元数乘以四元数的逆
四元数的加减
负四元数
四元数的模 Quat::length()
四元数的共轭 Quat::conj()
四元数的逆 Quat::inverse() = Quat::conj() / Quat::length2()
与旋转相关的操作:这部分是整个四元数最为关键的部分,其作用就体现于此,通过重载makeRotate()来实现四元数与矩阵和欧拉角的转换
四元数能被解释为角位移的轴――角对方式 q = (nxsin(△/2), nysin(△/2),nzsin(△/2), cos (△/2))
在OSG中,makeRotate()还实现是通过三个轴--角对来生成一个四元数的方法,是通过每个轴――角对生成一个四元数,然后叉乘三个四元数得到一个四元数。
makeRotate()还实现了获取将一个向量旋转到另一个向量的四元数,其关键是获取旋转角(向量点乘)和旋转轴(向量叉乘)。
将四元数转换为轴-角对
四元数插值—slerp:slerp是3D数学中四元数存在理由,它是球面线性插值的缩写(Spherical Linear Interpolotioin),由此可看出其重要性。Slerp运算可以在两个四元数之间平滑插值,其运算避免了欧拉角插值的所有问题。
用四元数实现旋转向量:这个方法采用了nVidia SDK的实现方法
以上是OSG::Quat中实现的四元数的方法,但四元数还有其它的一些运算,也值得一提,所以将实现为四元数的外部函数
四元数求幂:该操作可以从角位移中抽取“一部分” Quat::pow()
单位四元数:Quat::normalize()
3. 实际应用
我们可以扩展一个标准的3D点(x,y,z)到四元数空间,通过定义四元数p=[0,(x, y, z)]即可。
4. 疑惑
makeRotate(const Vec3d& from, const Vec3d& to)中的算法不胜理解,makeRotate_original()方法的目的与算法也不太明白
代码实例(http://blog.csdn.net/pizi0475/article/details/5584610)- // HPR to Quat
- osg::Quat HPRToQuat(double heading, double pitch, double roll)
- {
- osg::Quat q(
- roll,osg::Vec3d(0.0, 1.0, 0.0),
- pitch,osg::Vec3d(1.0, 0.0, 0.0),
- heading,osg::Vec3d(0.0, 0.0, 1.0));
- return q;
- }
- // Quat to HPR,pitch范围:[-PI/2, PI/2]
- void QuatToHPR(osg::Quat q, double& heading, double& pitch, double& roll)
- {
- double test = q.y() * q.z() + q.x() * q.w();
- if (test > 0.4999)
- { // singularity at north pole
- heading = 2.0 * atan2(q.y(), q.w());
- pitch = osg::PI_2;
- roll = 0.0;
- return;
- }
- if (test < -0.4999)
- { // singularity at south pole
- heading = 2.0 * atan2(q.y(), q.w());
- pitch = -osg::PI_2;
- roll = 0.0;
- return;
- }
- double sqx = q.x() * q.x();
- double sqy = q.y() * q.y();
- double sqz = q.z() * q.z();
- heading = atan2(2.0 * q.z() * q.w() - 2.0 * q.y() * q.x(), 1.0 - 2.0 * sqz - 2.0 * sqx);
- pitch = asin(2.0 * test);
- roll = atan2(2.0 * q.y() * q.w() - 2.0 * q.z() * q.x(), 1.0 - 2.0 * sqy - 2.0 * sqx);
- }
osg::Quat 学习心得,和暂时迷惑者共同探讨
osg:: Quat的心得:
先介绍四元数:Q = [w,(x,y,z)]被定义为一个四元数,w为一个实数,(x,y,z)是一个三维矢量,四元数的基底为(1,i,j,k),则 Q = w + xi + yj + zk;四元数是复数在四维空间的推广,于是可以认为i,j,k是四元数的虚单位。
欧拉证明了任何旋转都可以通过对一个轴的单一旋转来表示,因此方向可以分解为一个角度(θ)+轴(n),可能有人和我一样以为osg:: Quat(x,y,z,w)就是沿着轴(x,y,z)正向看去逆时针旋转w角度(正向是从原点指向点(x,y,z)的矢量的方向)但是,这是不对的。欧拉认为这个角度和轴不是简单的放在四元数中,他们的关系应该如下:
q = [ cos(θ/2) , sin(θ/2)n ] = [ cos(θ/2), ( sin(θ/2)x, sin(θ/2)y, sin(θ/2)z ) ]
四元数的标量部分为cos(θ/2),即把旋转角的一半的余弦值作为四元数的标量部分;把旋转角的一半的正弦值和这个旋转轴数乘以后得到的矢量作为四元数的矢量部分。而osg:: Quat在处理时将标量部分放在了最后,上面q用在osg中用osg:: Quat表示为:
q = [ sin(θ/2)n , cos(θ/2) ] = [ ( sin(θ/2)x, sin(θ/2)y, sin(θ/2)z ) ,cos(θ/2) ]
即在osg:: Quat中,它的成员变量_v[3]表示的是四元数的标量部分,因为余弦函数的值域是有界闭区间[-1,1],所以看来_v[3]的值是介于[-1,1]区间上的,看osg:: Quat的一个成员函数makeRotate:
void Quat::makeRotate( value_type angle, value_type x, value_type y, value_type z )
{
const value_type epsilon = 0.0000001;
value_type length = sqrt( x*x + y*y + z*z );
if (length < epsilon)
{
// ~zero length axis, so reset rotation to zero.
*this = Quat();
return;
}
value_type inversenorm = 1.0/length; //等会用于单位化旋转轴
value_type coshalfangle = cos( 0.5*angle ); //旋转角一半的余弦,用于实部
value_type sinhalfangle = sin( 0.5*angle ); //旋转角一半的正弦
_v[0] = x * sinhalfangle * inversenorm;
_v[1] = y * sinhalfangle * inversenorm;
_v[2] = z * sinhalfangle * inversenorm; //旋转轴被单位化
_v[3] = coshalfangle;
}
函数最后一句 "_v[3] = coshalfangle;" 说明了osg:: Quat对实部处理的这个事实,而同时,旋转轴被单位化,再乘以一个值域有界为[-1,1]的正弦值,所以虚部每一个分量都介于[-1,1]之间。通过上述处理,可以将一个旋转角,一个旋转轴转化为一个四元数。
再看看getRotate这个函数,怎样从一个四元数得到它的旋转角度和转轴?
void Quat::getRotate( value_type& angle, value_type& x, value_type& y, value_type& z ) const
{
value_type sinhalfangle = sqrt( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] );
angle = 2.0 * atan2( sinhalfangle, _v[3] );
if(sinhalfangle)
{
x = _v[0] / sinhalfangle;
y = _v[1] / sinhalfangle;
z = _v[2] / sinhalfangle;
}
else
{
x = 0.0;
y = 0.0;
z = 1.0;
}
}
代码中第一句 "sqrt( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] )" 为什么会是旋转角一半的正弦值,看看刚才makeRotate的时候是怎么处理过的:
_v[0] = x * sinhalfangle * inversenorm;
_v[1] = y * sinhalfangle * inversenorm;
_v[2] = z * sinhalfangle * inversenorm;
代入上式,提取公因子sinhalfangle,可以发现正确。因为向量部分已经单位化了,所以可以直接得到时sin(angle/2)。事实上,_v[3]就是cos(angle/2),所以也可以
sin(angle/2) = sqrt(1 - _v[3]*_v[3]),_v[3]就是cos(angle/2),所以angle = 2 * acos(_v[3]),osg:: Quat中旋转角angle的得到比较迂回,它求得coshalfangle/sinhalfangle的反正切,即得angle的一半。求旋转轴的时候分别给_v[0],_v[1],_v[2]除以sinhalfangle得到单位化后的旋转轴。
四元数的运算法则:i*i = j*j = k*k = i*j*k = -1;
ij = k; ji = -k;
jk = i; kj = -i;
ki = j; ik = -j;
大家可以用上面的法则推倒一下osg:: Quat中四元数的乘法,发现结果是正确的^^,乘法的时候,osg:: Quat是将实部放在最后一个分量_v[3]上的,有四元数的共轭:
inline Quat conj () const { return Quat( -_v[0], -_v[1], -_v[2], _v[3] ); }可以进一步看出第四个分量是实部,欢迎讨论,批评指正。