四元数入门

4元数宝典

这是国内找不到的超好文章。（为什么大陆的4元数文章很垃圾呢？）
（翻译中。。。奉献给大家~~）

70秒即懂，能使用，用四元数，4元数，阔特尼恩，Quaternion旋转
(C) 中田  亨  （独立行政法人  产业技术综合研究所  数字人类研究中心  研究员 博士（工学）） 
2003年11月25日

★这个页面的对象读者
想把三次元的旋转，用CG等定量地处理的人
使用欧拉角(Euler Angles)的话，不懂得其道理的人 
卡尔丹角和欧拉角(Cardan Angles)不能区别的人 
对吉恩瓦尔洛克很困惑的人
但是，对数学之类麻烦的事情很讨厌的人
想要实例程序的人 
没有时间的人

★旋转篇：
　我将说明使用了四元数（si yuan shu, quaternion）的旋转的操作步骤
（１）四元数的虚部，实部和写法
所谓四元数，就是把4个实数组合起来的东西。
4个元素中，一个是实部，其余3个是虚部。
比如，叫做Q的四元数，实部t而虚部是x,y,z构成，则像下面这样写。
Q = (t; x, y, z) 
又，使用向量 V=(x,y,z)，
Q = (t; V)  
也可以这么写。

正规地用虚数单位i,j,k的写法的话，
Q = t + xi + yj + zk 
也这样写，不过，我不大使用

（２）四元数之间的乘法
虚数单位之间的乘法 
ii = -1, ij = -ji = k (其他的组合也是循环地以下同文) 
有这么一种规则。（我总觉得，这就像是向量积（外积），对吧） 
用这个规则一点点地计算很麻烦，所以请用像下面这样的公式计算。

A = (a; U) 
B = (b; V) 
AB = (ab - U·V; aV + bU + U×V)
不过，“U·V”是内积，「U×V」是外积的意思。
注意：一般AB<>BA所以乘法的左右要注意！

（3）3次元的坐标的四元数表示
如要将某坐标(x,y,z)用四元数表示，
P = (0; x, y, z) 
则要这么写。
 
另外，即使实部是零以外的值，下文的结果也一样。用零的话省事所以我推荐。

（４）旋转的四元数表示
以原点为旋转中心，旋转的轴是(α, β, γ)
（但 α^2 + β^2 + γ^2 = 1）， 
（右手系的坐标定义的话，望向向量(α, β, γ)的前进方向反时针地） 
转θ角的旋转，用四元数表示就是，
Q = (cos(θ/2); α sin(θ/2), β sin(θ/2), γ sin(θ/2)) 
R = (cos(θ/2); -α sin(θ/2), -β sin(θ/2), -γ sin(θ/2)) 
(另外R 叫 Q 的共轭四元数。）

那么，如要实行旋转，
则 R P Q = (0; 答案)

请像这样三明治式地计算。这个值的虚部就是旋转之后的点的坐标值。
 （另外，实部应该为零。请验算看看）

*未完。。。

instemast_REAL 2007-2-24 19:04:57
 
 注册: 2007-2    状态: Offline 1   Top 
 
instemast_REAL
Exp:124

侦察兵
  发表于: 2007-2-24 19:06:00 档案 | 短信 | 树状 | 收藏 | 编辑 | 删除 | 引用   
 

--------------------------------------------------------------------------------
Re:4元数宝典

/// Quaternion.cpp 
/// (C) Toru Nakata, toru-nakata@aist.go.jp 
/// 2004 Dec 29 
  
#include <math.h> 
#include <iostream.h> 
  
/// Define Data type 
typedef struct 
{ 
              double t; // real-component 
              double x; // x-component 
              double y; // y-component 
              double z; // z-component 
} quaternion; 
  
  
Kakezan 
quaternion Kakezan(quaternion left, quaternion right) 
{ 
              quaternion ans; 
              double d1, d2, d3, d4; 
  
              d1 =  left.t * right.t; 
              d2 = -left.x * right.x; 
              d3 = -left.y * right.y; 
              d4 = -left.z * right.z; 
              ans.t = d1+ d2+ d3+ d4; 
  
              d1 =  left.t * right.x; 
              d2 =  right.t * left.x; 
              d3 =  left.y * right.z; 
              d4 = -left.z * right.y; 
              ans.x =  d1+ d2+ d3+ d4; 
  
              d1 =  left.t * right.y; 
              d2 =  right.t * left.y; 
              d3 =  left.z * right.x; 
              d4 = -left.x * right.z; 
              ans.y =  d1+ d2+ d3+ d4; 
  
              d1 =  left.t * right.z; 
              d2 =  right.t * left.z; 
              d3 =  left.x * right.y; 
              d4 = -left.y * right.x; 
              ans.z =  d1+ d2+ d3+ d4; 
              
              return ans; 
} 
  
Make Rotational quaternion 
quaternion MakeRotationalQuaternion(double radian, double AxisX, double AxisY, double AxisZ) 
{ 
              quaternion ans; 
              double norm; 
              double ccc, sss; 
              
              ans.t = ans.x = ans.y = ans.z = 0.0; 
  
              norm = AxisX *  AxisX +  AxisY *  AxisY +  AxisZ *  AxisZ; 
              if(norm <= 0.0) return ans; 
  
              norm = 1.0 / sqrt(norm); 
              AxisX *= norm; 
              AxisY *= norm; 
              AxisZ *= norm; 
  
              ccc = cos(0.5 * radian); 
              sss = sin(0.5 * radian); 
  
              ans.t = ccc; 
              ans.x = sss * AxisX; 
              ans.y = sss * AxisY; 
              ans.z = sss * AxisZ; 
  
              return ans; 
} 
  
Put XYZ into  quaternion 
quaternion PutXYZToQuaternion(double PosX, double PosY, double PosZ) 
{ 
              quaternion ans; 
  
              ans.t = 0.0; 
              ans.x = PosX; 
              ans.y = PosY; 
              ans.z = PosZ; 
  
              return ans; 
} 
  
/ main 
int main() 
{ 
              double px, py, pz; 
              double ax, ay, az, th; 
              quaternion ppp, qqq, rrr; 
  
              cout << "Point Position (x, y, z) " << endl; 
              cout << "  x = "; 
              cin >> px; 
              cout << "  y = "; 
              cin >> py; 
              cout << "  z = "; 
              cin >> pz; 
              ppp = PutXYZToQuaternion(px, py, pz); 
  
              while(1) { 
                            cout << "/nRotation Degree ? (Enter 0 to Quit) " << endl; 
                            cout << "  angle = "; 
                            cin >> th; 
                            if(th == 0.0) break; 
  
                            cout << "Rotation Axis Direction ? (x, y, z) " << endl; 
                            cout << "  x = "; 
                            cin >> ax; 
                            cout << "  y = "; 
                            cin >> ay; 
                            cout << "  z = "; 
                            cin >> az; 
  
  
                            th *= 3.1415926535897932384626433832795 / 180.0; /// Degree -> radian; 
  
                            qqq = MakeRotationalQuaternion(th, ax, ay, az); 
                            rrr = MakeRotationalQuaternion(-th, ax, ay, az); 
  
                            ppp = Kakezan(rrr, ppp); 
                            ppp = Kakezan(ppp, qqq); 
  
                            cout << "/nAnser X = " << ppp.x 
                                          <<  "/n      Y = " << ppp.y 
                                          <<  "/n      Z = " << ppp.z << endl; 
  
              } 
  
              return 0; 
} 

*未完。。。

 
 

四元数的长处和缺点

冰川球的吉恩瓦尔。

吸收前后左右的倾斜而保持水平

长处

“直观的”

如果给出旋转轴和旋转角度的话，就算不考虑欧拉角之类的，也可以立即计算。

“连续性”

表示相似的旋转的四元数的值也相似。没有吉恩瓦尔洛克现象。
（所谓吉恩瓦尔洛克就是，“北极和南极等特殊点，在自转的情况下是不能动的。”“在中国餐馆的圆桌上，
放在正中央的酱油怎么转也转不近”等，有效的自由度和次元丧失的现象。在软件中成为例外处理的困难因而是麻烦事。本来吉恩瓦尔是把罗盘针吊起来的机构。为了使得即便船摇晃罗盘针也保持水平而不摇晃。吉恩瓦尔的柄把罗盘针晃晃荡荡地吊着，但若使得这些东西的方向一致了，就不能晃了（=洛克lock），振动传到罗盘针了。）

“记忆效率好，计算快”

旋转能够只用4个数值记述。与阿弗因变换矩阵等相比，不用记忆。减少不必要的计算。
 
缺点

“从外观上看不知道表示什么意思”
 
从四元数的成分来看，其表示什么意思，一眼看不出来。
 
“不能表示多圈旋转”
如你所见，使用cos和sin，所以θ是10度，370度，还是-350度，不能区别。
想要制作咕噜咕噜转许多圈的动画的情况下，不想只旋转一遍，请仔细区分。

“是以原点为中心的旋转”
想以原点以外为中心的情况，请让坐标穿上木屐，旋转，然后还原木屐。
（下文的一般位移篇也check一下吧！）

四元数旋转的软件的源代码
<<source code 贴过了>>

一般位移篇：作为旋转以外的变形的扩大缩小（与并进）的步骤

考虑四元数的标量倍，就可以将旋转篇的结果坐标给扩大缩小。比如，
P = (0; x, y, z)

（另外实部可以为任意值。0好算）

kQ = (kcos(θ/2); kα sin(θ/2), kβ sin(θ/2), kγ sin(θ/2))

kR = (kcos(θ/2); -kα sin(θ/2), -kβ sin(θ/2), -kγ sin(θ/2))
　 
（KR是KQ的共轭四元数。）

则有，

kR P kQ = ( 0; 旋转后的XYZ坐标的k的自乘倍)

这样，就可以将以原点为中心的旋转，和以原点为中心的扩大缩小一起操作了。
这表示，可以自由地操作以原点为准的方位及其距离。就是说，将一点(X,Y,Z)映射到任意的点(X',Y',Z')的变换，可以用四元数表示。但X=Y=Z=0的情况例外。

制作你所希望的四元数篇：实施某旋转的四元数的作法

表示将一点A，移到别的点B的以原点为中心的旋转的四元数，像这样制作。

旋转轴（α，β，γ）和位置向量的外积 Ｂ×Ａ 平行且同向。（A和B的顺序，按四元数R和Q的乘法的顺序交替。本文的乘法顺序的定义的话，就是BA的顺序。）作了除法之后，将把长度改为1。

（α，β，γ）＝ (Ｂ×Ａ) / (｜Ｂ×Ａ｜)

cos(θ)是，将位置向量的内积 Ｂ·Ａ 除以｜Ｂ｜·｜Ａ｜。
cos(θ)＝ （Ｂ·Ａ） / (｜Ｂ｜·｜Ａ｜)

根据三角函数的半角公式，
cos(θ/2)＝±√{ 0.5 * (1 + cosθ) } 
sin(θ/2)＝±√{ 0.5 * (1 - cosθ) }

正负号呢，由于 0°≦θ≦180°所以，cos(θ/2)≧０ 且 sin(θ/2)≧０ 就是说，两边都为正。
cos(θ/2)＝√{ 0.5 * (1 + cosθ) } 
sin(θ/2)＝√{ 0.5 * (1 - cosθ) }

也要进行扩大缩小的情况、k = √{｜Ｂ｜/｜Ａ｜}

从欧拉角到四元数的变换及其逆变换的问题（背景说明篇）

这个问题时不时被谈论，由于数学原因和数学之外的事情而很难。
 
3次元的旋转呢，
(1) 作为旋转轴的直线
(2) 旋转角度
有这两个信息，就可以表示。这便是“万能的3次元旋转表示”。

把这样的单次的旋转叫做single rotation。

（注意：single rotation 没有旋转中心点。只有旋转轴。如果是，关于像球形关节这样的旋转的，拥有不同旋转轴的single rotation 复合而成的 rotation sequence ，的话，就能够考虑旋转中心点了。）
（*译注：此句的复句结构复杂，但无论怎么理解都不影响意思！）
 
（注意：四元数将single ratation以旋转轴和旋转角度的信息来直接地表示。
欧拉角呢，是要把一个single rotation，用分解为3个single rotation的方式来表示。可以说故意搞复杂了。而且，蕴含着引起吉恩瓦尔洛克的问题的可能性。）

“欧拉角-->四元数”
考虑像这样的直接的变换是非常难的问题。
正确的是，应该取，
像“欧拉角-->万能表示-->四元数”
这样的，
“一种表式形式-->万能表示-->别的表示形式”
的路径。我想这最安全也易思考。

 

文章主要内容在前面，后面闲谈比较多。而且后面的语法有点乱（科学家语文不好，还喜欢用复杂句式）
但注意，他告诉了我们什么是真正的欧拉角，这有点出乎意料，和一般d3d书上说的不一样。

变换困难的数学原因

四元数：原点不在旋转轴上的情况，必须将其木屐（偏移量）用别的途径表示。
欧拉角：不仅有原点的偏移量的问题，而且，产生吉恩瓦尔洛克现象情况下必须做例外处理。还有，轴的旋转顺序（虽说有惯例但）不统一。（总的来说欧拉角有太多使用上的难点。）

变换困难的数学之外的原因

数据格式因软件而异。“一段变换程序”只在一个公司的软件范围内通用。 
坐标轴的设定不同。右手系还是左手系不统一。还有“向上”和“向右”是哪个轴等问题。也有历史原因，所以不统一。
角度的正负方向的决定方式难以知晓。数据格式与坐标轴设定的技术，光看源代码不知道。
(特别是，拥有三个角度的欧拉角很麻烦)
既如此，为什么要用欧拉角呢？
because there are many gyros whitch are roll,pitch and yaw,in the output?
像大炮的底座一样地，机构使用了roll,pitch,yaw的情况也是有的。
(注：gyro = gyrocompass旋转罗盘,gyroscope陀螺仪)

欧拉角是什么

Euler考虑到，“一个三次元的single rotation，可以以组合三个以坐标轴作为旋转轴的三次元旋转的方式来表示”。这是广义的欧拉角。
就是说，以，从坐标系看来倾斜的轴，为旋转轴的旋转，搞得好的话，首先以Z轴中心旋转α度，然后以X轴中心旋转β度，最后以Z轴中心旋转γ度，以这种方式，是可以作出一模一样的旋转的。（发生吉恩瓦尔洛克现象的情况，失去了解的一义性，成为了件麻烦的事情。）
欧拉自己，把旋转轴的选择方式使用了Z-X-Z。这个选轴顺序的情况，是狭义的欧拉角。（考虑陀螺的运动的时候，首先把陀螺的高速自转用Z轴旋转表示，把陀螺的芯棒倒向地面的旋转用X轴旋转表示把陀螺芯棒的水平旋转用Z轴旋转表示。）

世上像Z-Y-Z、X-Z-X, Y-X-Y之类的模式，似乎也有人称之为欧拉角的。
轴的选择方式，也有，像X-Y-Z，Y-Z-Z这样的3轴顺次使用的模式。（另外，3个轴的立场不是平等的。存在着由旋转操作的先后所引起的不同

。）也有人把这些认作是广义的欧拉角，而，只承认狭义的欧拉角的人是误解之源。

X-Y-Z式是被频繁使用的模式。这被叫做roll,pitch,yaw角式，也有人叫卡尔丹角(Cardan angles)。

从欧拉角到四元数表示

这是极其简单的，把欧拉角的3次旋转，分别用四元数表示，顺次作乘法即可。

从四元数到欧拉角

「オイラー角→四元数」の逆変換を考えることになります。公式化してもいいのですが、オイラー角がZ-X-Z式なのかX-Y-Z式なのかの違い

で、大きく変わってしまいます。シンバルロック発生時には、公式の中でゼロ割りが起きますので、その処理も面倒みてください。
シングルローテーションの不思議
我们来考虑一下“欧拉角-〉四元数”的逆变换。虽说可以公式化但是，因欧拉角是Z-X-Z式还是X-Y-Z式的不同，就会有大的变化。吉恩瓦尔洛

克现象发生的时候，在公式中会发生零作除数的情况，其处理也麻烦请看。
single rotation的不可思议
一个single rotation，可以用一个四元数来表示的。
于是可以说，进行N次single rotation的旋转，能够用N个四元数的积的一个四元数来表示。也就是说，多个single rotation的串珠实行，无

视其过程，只看其结果的话，可以只用1个四元数表示。
结果，就算咕噜咕噜地转，也只不过是进行了1个single rotation。就算把足球转到那边，再转回这边，按照初期姿势和最终姿势，完全没有

移动的点至少有2个！或者说，在1次旋转中，到最终姿势，只是跳个一下就可以了。