OpenGL基本瞭解(九) 四元數在OPENGL的應用

Quaternion 四元數在OPENGL的應用

依照 Wiki 上的一篇介紹 Rotation/Orientation 表示法的 文章說法，要描述一個物體在三維空間內的旋轉(Rotation)/面向(Orientation)可以用很多不同的代數系統來表示。其中在 3D 繪圖裡最常用的三種是矩陣系統(Matrix)、尤拉系統(Euler Axis & Angles)、以及四元數系統(Quaternion)。

尤拉系統用於旋轉時會遇到一個很不幸的問題 - 萬向節鎖 (Gimbal Lock)。一篇說明什麼是萬向節鎖的文章

三維動畫師最厭惡的情況之一『萬向鎖（Gimbal Lock）』就是這個傢伙的問題。按不同軸以優拉角旋轉幾次後，出現x,y,z三個軸完全變成同向的情況，也就是說，優拉角很容易出現旋轉到最後只剩一個方向可以旋轉的情況，這就是恐怖的『萬向鎖』。

有一個影片可以作為說明  https://www.youtube.com/watch?v=rsKy-4dbA04

要避免萬向節鎖，一個方法為是轉換成四元數系統。


四元數是一個很深奧的數學概念，用來解決尤拉系統的旋轉問題只是它廣大應用中的一小部分。所以若直接去查閱有關四元數的教學文章是很吃力且發散的。兩篇文章不錯的說明: 一篇是 Wiki 上專門講解四元數用於旋轉問題的文章，另一篇是 GameDev 上介紹四元數的實例用途的文章。 建議以 GameDev 那篇為主、Wiki 那篇為輔去閱讀，因為 Wiki 那篇還是有一些數學導證，對於程式開發人員來說比較不是那麼重要。GameDev 那篇裡面提到的那些不同系統間的轉換公式，都是以 OpenGL 使用的右手座標系統為準去推導的。Wiki 那篇後面有一個章節比較了三個表示法在不同操作時的效能差異。很有趣的是，在進行連續多組旋轉互相乘積(組合)時，四元數所需的運算量比矩陣少的; 但是在進行與點乘積(取變換後位置)時，因為四元數還是要先轉換回矩陣才有辦法乘，所以是矩陣佔優勢。

有關四元數的資料可參考 wiki >>  http://zh.wikipedia.org/wiki/%E5%9B%9B%E5%85%83%E6%95%B8

在3D運算中的使用可參考 http://blog.csdn.net/kesalin/article/details/2187347


但四元數不是『絕對完美』，因為插值的時候過渡速率不恆定，且很難解決。不過這比起『恐怖萬向鎖』已經是很小的問題。


關於參考的OPENGL數學運算代碼可參考 
http://www.linuxgraphics.cn/opengl/opengl_quaternion.html

http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/index.htm


對於3D的旋轉可參考Yaw, Pitch, Roll的含义 

主要參考網頁
1.   http://blog.roodo.com/sayaku/archives/19544672.html

2.   http://coazure-code.blogspot.tw/2010/01/pv3d.html

3.   http://www.wretch.cc/blog/zevoid/658309



以下是參考程式碼的轉貼

#pragma once

#include "Vector.h"
#include "GLESMath.h"

struct Quaternion
{
    float x;
    float y;
    float z;
    float w;
   
    Quaternion();
    Quaternion(float x, float y, float z, float w);
   
    Quaternion Slerp(float mu, const Quaternion& q) const;
    Quaternion Rotated(const Quaternion& b) const;
    Quaternion Scaled(float scale) const;
   
    float Dot(const Quaternion& q) const;
    void ToMatrix4(KSMatrix4 * m) const;
    Vector4<float> ToVector() const;
    void ToIdentity();
   
    Quaternion operator-(const Quaternion& q) const;
    Quaternion operator+(const Quaternion& q) const;
    bool operator==(const Quaternion& q) const;
    bool operator!=(const Quaternion& q) const;
   
    void Normalize();
    void Rotate(const Quaternion& q);
   
    static Quaternion CreateFromVectors(const Vector3<float>& v0, const Vector3<float>& v1);
    static Quaternion CreateFromAxisAngle(const Vector3<float>& axis, float radians);
};

inline Quaternion::Quaternion() : x(0), y(0), z(0), w(1)
{}

inline Quaternion::Quaternion(float x, float y, float z, float w) : x(x), y(y), z(z), w(w)
{}

inline void Quaternion::ToIdentity()
{
    x = y = z = 0;
    w = 1.0;
}

// Ken Shoemake's famous method.
inline Quaternion Quaternion::Slerp(float t, const Quaternion& v1) const
{
    const float epsilon = 0.0005f;
    float dot = Dot(v1);
   
    if (dot > 1 - epsilon) {
        Quaternion result = v1 + (*this - v1).Scaled(t);
        result.Normalize();
        return result;
    }
   
    if (dot < 0)
        dot = 0;
   
    if (dot > 1)
        dot = 1;
   
    float theta0 = acos(dot);
    float theta = theta0 * t;
   
    Quaternion v2 = (v1 - Scaled(dot));
    v2.Normalize();
   
    Quaternion q = Scaled(cos(theta)) + v2.Scaled(sin(theta));
    q.Normalize();
    return q;
}

inline Quaternion Quaternion::Rotated(const Quaternion& b) const
{
    Quaternion q;
    q.w = w * b.w - x * b.x - y * b.y - z * b.z;
    q.x = w * b.x + x * b.w + y * b.z - z * b.y;
    q.y = w * b.y + y * b.w + z * b.x - x * b.z;
    q.z = w * b.z + z * b.w + x * b.y - y * b.x;
    q.Normalize();
    return q;
}

inline Quaternion Quaternion::Scaled(float s) const
{
    return Quaternion(x * s, y * s, z * s, w * s);
}

inline float Quaternion::Dot(const Quaternion& q) const
{
    return x * q.x + y * q.y + z * q.z + w * q.w;
}

inline void Quaternion::ToMatrix4(KSMatrix4 * result) const
{
    const float s = 2;
    float xs, ys, zs;
    float wx, wy, wz;
    float xx, xy, xz;
    float yy, yz, zz;
    xs = x * s;  ys = y * s;  zs = z * s;
    wx = w * xs; wy = w * ys; wz = w * zs;
    xx = x * xs; xy = x * ys; xz = x * zs;
    yy = y * ys; yz = y * zs; zz = z * zs;
   
    result->m[0][0] = 1 - (yy + zz);
    result->m[0][1] = xy + wz;
    result->m[0][2] = xz - wy;
    result->m[0][3] = 0;
   
    result->m[1][0] = xy - wz;
    result->m[1][1] = 1 - (xx + zz);
    result->m[1][2] = yz + wx;
    result->m[1][3] = 0;
   
    result->m[2][0] = xz + wy;
    result->m[2][1] = yz - wx;
    result->m[2][2]= 1 - (xx + yy);
    result->m[2][3] = 0;
   
    result->m[3][0] = 0;
    result->m[3][1] = 0;
    result->m[3][2] = 0;
    result->m[3][3] = 1;
}

inline Vector4<float> Quaternion::ToVector() const
{
    return Vector4<float>(x, y, z, w);
}

inline Quaternion Quaternion::operator-(const Quaternion& q) const
{
    return Quaternion(x - q.x, y - q.y, z - q.z, w - q.w);
}

inline Quaternion Quaternion::operator+(const Quaternion& q) const
{
    return Quaternion(x + q.x, y + q.y, z + q.z, w + q.w);
}

inline bool Quaternion::operator==(const Quaternion& q) const
{
    return x == q.x && y == q.y && z == q.z && w == q.w;
}

inline bool Quaternion::operator!=(const Quaternion& q) const
{
    return !(*this == q);
}

inline void Quaternion::Normalize()
{
    *this = Scaled(1 / sqrt(Dot(*this)));
}

inline void Quaternion::Rotate(const Quaternion& q2)
{
    Quaternion q;
    Quaternion& q1 = *this;
   
    q.w = q1.w * q2.w - q1.x * q2.x - q1.y * q2.y - q1.z * q2.z;
    q.x = q1.w * q2.x + q1.x * q2.w + q1.y * q2.z - q1.z * q2.y;
    q.y = q1.w * q2.y + q1.y * q2.w + q1.z * q2.x - q1.x * q2.z;
    q.z = q1.w * q2.z + q1.z * q2.w + q1.x * q2.y - q1.y * q2.x;
   
    q.Normalize();
    *this = q;
}

// Compute the quaternion that rotates from a to b, avoiding numerical instability.
// Taken from "The Shortest Arc Quaternion" by Stan Melax in "Game Programming Gems".
//
inline Quaternion Quaternion::CreateFromVectors(const Vector3<float>& v0, const Vector3<float>& v1)
{
    if (v0 == -v1)
        return Quaternion::CreateFromAxisAngle(vec3(1, 0, 0), Pi);
   
    Vector3<float> c = v0.Cross(v1);
    float d = v0.Dot(v1);
    float s = sqrt((1 + d) * 2);
   
    Quaternion q;
    q.x = c.x / s;
    q.y = c.y / s;
    q.z = c.z / s;
    q.w = s / 2.0f;
   
    return q;
}

inline Quaternion Quaternion::CreateFromAxisAngle(const Vector3<float>& axis, float radians)
{
    Quaternion q;
    q.w = cos(radians / 2);
    q.x = q.y = q.z = sin(radians / 2);
    q.x *= axis.x;
    q.y *= axis.y;
    q.z *= axis.z;
   
    return q;
}