我将说明使用了四元数(si yuan shu, quaternion)的旋转的操作步骤
( 1 )四元数的虚部,实部和写法
所谓四元数,就是把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
也这样写,不过,我不大使用
( 2 )四元数之间的乘法
虚数单位之间的乘法
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)
则要这么写。
另外,即使实部是零以外的值,下文的结果也一样。用零的话省事所以我推荐。
( 4 )旋转的四元数表示
以原点为旋转中心,旋转的轴是(α, β, γ)
(但 α ^ 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;
}
转自:
http://www.cppblog.com/tommy/archive/2010/10/21/130772.html?opt=admin