Posted on 2010-10-21 17:26
Tommy Liang 閱讀(8560)
評論(0) 編輯 收藏 引用
四元數常常可以在3D的書上看到。
但我的那本3D圖形學書上，在沒講四元數是干什么的之前，就列了幾張紙的公式，
大概因為自己還在上高中，不知道的太多，看了半天沒看懂。。。
終于，在gameres上看到了某強人翻譯的一個“4元數寶典 ”（原文是日本人寫的。。。），感覺很好，分享下。
★旋轉篇：
我將說明使用了四元數（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; 答案)
請像這樣三明治式地計算。這個值的虛部就是旋轉之后的點的坐標值。
（另外，實部應該為零。請驗算看看）
例子代碼
/// 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;
//// Bill 注：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;
}