深入淺出四元數（誤）

作者：~XD│2018-03-31 22:57:11│巴幣：4│人氣：569

看到有人問這個，就稍微了解一下下，順便留作記錄…

原串：https://forum.gamer.com.tw/C.php?bsn=60602&snA=2054

不過昨天看了你的文還是跑去了解一下四元數
找到這個網頁是我個人覺得幫助最大的教學
http://www015.upp.so-net.ne.jp/notgeld/quaternion.html
（別告北為什麼是日文，找不到好的中文教學我也很絕望）

裡面有段四元數的程式碼，看了大概就知道四元數是在搞什麼

/// 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;
}

我試跑程式得到下面那樣，有多加一些cout方便觀察東西是怎麼變出來的

Point Position (x, y, z)
  x = 1
  y = 0
  z = 0

Rotation Degree ? (Enter 0 to Quit)
  angle = 90
Rotation Axis Direction ? (x, y, z)
  x = 0
  y = 0
  z = 1

qqq   T = 0.707107
      X = 0
      Y = 0
      Z = 0.707107

rrr   T = 0.707107
      X = -0
      Y = -0
      Z = -0.707107

r*p   T = 0
      X = 0.707107
      Y = -0.707107
      Z = 0

r*p*q T = 0
      X = 2.22045e-016
      Y = -1
      Z = 0

Rotation Degree ? (Enter 0 to Quit)
  angle =

四元數結構(t,x,y,z)
我拿點座標(1,0,0)當作參數，直接複製成四元數(0,1,0,0)
接著旋轉90度，旋轉軸我取(0,0,1)
先做單位化然後應用這個公式
Q = (cos(θ/2); α sin(θ/2), β sin(θ/2), γ sin(θ/2))
R = (cos(θ/2); -α sin(θ/2), -β sin(θ/2), -γ sin(θ/2))

角度取半，r取負角，q取正角
實部t取cos
虛部xyz取sin乘上對應的來源坐標分量
本來這段公式打成方程式一看就明，列這樣真是難看

得到q和r
最後算r*p*q

這個乘法是四元數乘法，計算規則是
A = (a; U)
B = (b; V)
AB = (ab - U・V; aV + bU + U×V)
內有一個點積一個叉積
實際計算你自己看quaternion Kakezan(quaternion left, quaternion right)就明白了

然後得到(0,x,y,z)，裡面的xyz就是旋轉後的新座標
完成四元數運算

做一次四元數計算大概需要
產生一對共軛四元數rq，有兩次單位化和sin、cos計算
p為要被旋轉的點座標參數，計算r*p*q，有兩次四元數乘法計算

四元數對歐拉角的優勢在我看來大概是如果要轉複雜的角度
歐拉角要對xyz軸轉三次，四元數可以一次轉三軸
四元數計算比純矩陣簡單，不然望向下面那個矩陣…

維基百科上面還有四元數的矩陣表達形式