最近要深入使用模型變換的東西，其中有個非常重要的就是旋轉，旋轉用經常用四元數來做，但是這個原理一直沒有弄懂，今天看資料，總結下：

（一）四元數:

   看原始碼中可以瞭解，網上別人總結的很好：

四元數：Q = [w,(x,y,z)]被定義為一個四元數,w為一個實數，(x,y,z)是一個三維向量，四元數的基底為(1,i,j,k),則 Q = w + xi + yj + zk;四元數是複數在四維空間的推廣，於是可以認為i,j,k是四元數的虛單位。
尤拉證明了任何旋轉都可以通過對一個軸的單一旋轉來表示，因此方向可以分解為一個角度(θ)+軸(n)，可能有人和我一樣以為osg:: Quat(x,y,z,w)就是沿著軸(x,y,z)正向看去逆時針旋轉w角度(正向是從原點指向點(x,y,z)的向量的方向)但是，這是不對的。尤拉認為這個角度和軸不是簡單的放在四元數中

void Quat::makeRotate( value_type angle, value_type x, value_type y, value_type z )
{
    const value_type epsilon = 0.0000001;

    value_type length = sqrt( x*x + y*y + z*z );
    if (length < epsilon)
    {
        // ~zero length axis, so reset rotation to zero.
        *this = Quat();
        return;
    }

    value_type inversenorm  = 1.0/length;            //等會用於單位化旋轉軸
    value_type coshalfangle = cos( 0.5*angle );      //旋轉角一半的餘弦，用於實部
    value_type sinhalfangle = sin( 0.5*angle );      //旋轉角一半的正弦

    _v[0] = x * sinhalfangle * inversenorm;
    _v[1] = y * sinhalfangle * inversenorm;
    _v[2] = z * sinhalfangle * inversenorm;          //旋轉軸被單位化
    _v[3] = coshalfangle;
}
 

void Quat::getRotate( value_type& angle, value_type& x, value_type& y, value_type& z ) const
{
    value_type sinhalfangle = sqrt( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] );

    angle = 2.0 * atan2( sinhalfangle, _v[3] );
    if(sinhalfangle)
    {
        x = _v[0] / sinhalfangle;
        y = _v[1] / sinhalfangle;
        z = _v[2] / sinhalfangle;
    }
    else
    {
        x = 0.0;
        y = 0.0;
        z = 1.0;
    }

}

程式碼中第一句 "sqrt( _v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] )" 為什麼會是旋轉角一半的正弦值，看看剛才makeRotate的時候是怎麼處理過的：
 _v[0] = x * sinhalfangle * inversenorm;
 _v[1] = y * sinhalfangle * inversenorm;
 _v[2] = z * sinhalfangle * inversenorm;

代入上式，提取公因子sinhalfangle，可以發現正確。因為向量部分已經單位化了，所以可以直接得到時sin(angle/2)。事實上，_v[3]就是cos(angle/2),所以也可以
sin(angle/2) = sqrt(1 - _v[3]*_v[3]),_v[3]就是cos(angle/2),所以angle = 2 * acos(_v[3]),osg:: Quat中旋轉角angle的得到比較迂迴，它求得coshalfangle/sinhalfangle的反正切，即得angle的一半。求旋轉軸的時候分別給_v[0],_v[1],_v[2]除以sinhalfangle得到單位化後的旋轉軸。



四元數的運演算法則：i*i = j*j = k*k = i*j*k = -1;
ij = k; ji = -k;
jk = i; kj = -i;
ki = j; ik = -j;
大家可以用上面的法則推倒一下osg:: Quat中四元數的乘法，發現結果是正確的^^,乘法的時候，osg:: Quat是將實部放在最後一個分量_v[3]上的，有四元數的共軛：
inline Quat conj () const { return Quat( -_v[0], -_v[1], -_v[2], _v[3] ); }可以進一步看出第四個分量是實部，歡迎討論，批評指正。

另一個函式：

void Quat::makeRotate( const Vec3d& from, const Vec3d& to )
{
    Vec3d sourceVector = from;
    Vec3d targetVector = to;


    value_type fromLen2 = from.length2();
    value_type fromLen;
    // normalize only when necessary, epsilon test
    if ((fromLen2 < 1.0-1e-7) || (fromLen2 > 1.0+1e-7)) {
        fromLen = sqrt(fromLen2);
        sourceVector /= fromLen;
    } else fromLen = 1.0;


    value_type toLen2 = to.length2();
    // normalize only when necessary, epsilon test
    if ((toLen2 < 1.0-1e-7) || (toLen2 > 1.0+1e-7)) {
        value_type toLen;
        // re-use fromLen for case of mapping 2 vectors of the same length
        if ((toLen2 > fromLen2-1e-7) && (toLen2 < fromLen2+1e-7)) {
            toLen = fromLen;
        }
        else toLen = sqrt(toLen2);
        targetVector /= toLen;
    }

//現在讓我們進入真正的東西使用“點積+ 1”ε測試,因為它可以重複使用    double dotProdPlus1 = 1.0 + sourceVector * targetVector;


    // 檢查退化的全面轉變。用ε來檢測
    if (dotProdPlus1 < 1e-7) {


        // 得到一個正交向量給定向量與最大平面向量的座標。然後用四元數與π軸角
        // Trick is to realize one value at least is >0.6 for a normalized vector.至少一個值> 0.6
        if (fabs(sourceVector.x()) < 0.6) {
            const double norm = sqrt(1.0 - sourceVector.x() * sourceVector.x());
            _v[0] = 0.0;
            _v[1] = sourceVector.z() / norm;
            _v[2] = -sourceVector.y() / norm;
            _v[3] = 0.0;
        } else if (fabs(sourceVector.y()) < 0.6) {
            const double norm = sqrt(1.0 - sourceVector.y() * sourceVector.y());
            _v[0] = -sourceVector.z() / norm;
            _v[1] = 0.0;
            _v[2] = sourceVector.x() / norm;
            _v[3] = 0.0;
        } else {
            const double norm = sqrt(1.0 - sourceVector.z() * sourceVector.z());
            _v[0] = sourceVector.y() / norm;
            _v[1] = -sourceVector.x() / norm;
            _v[2] = 0.0;
            _v[3] = 0.0;
        }
    }

    else {
        //找到最短的角度將歸一化向量轉換成另一個四元數。當向量共線的公式仍然有效
        const double s = sqrt(0.5 * dotProdPlus1);
        const Vec3d tmp = sourceVector ^ targetVector / (2.0*s);

//等到對應的四元數
        _v[0] = tmp.x();
        _v[1] = tmp.y();
        _v[2] = tmp.z();
        _v[3] = s;
    }
}

（二）四元數在旋轉中使用

後面就是對四元數的應用，主要還是應用與矩陣變換中，OSG中，矩陣直接對其使用Matrix（quat）

void Matrix_implementation::setRotate(const Quat& q)
{
    double length2 = q.length2();//四元數長度的平方（_v[0]*_v[0] + _v[1]*_v[1] + _v[2]*_v[2] + _v[3]*_v[3]）不能太小
    if (fabs(length2) <= std::numeric_limits<double>::min())
    {
        _mat[0][0] = 0.0; _mat[1][0] = 0.0; _mat[2][0] = 0.0;
        _mat[0][1] = 0.0; _mat[1][1] = 0.0; _mat[2][1] = 0.0;
        _mat[0][2] = 0.0; _mat[1][2] = 0.0; _mat[2][2] = 0.0;
    }
    else
    {
        double rlength2;
        // 如果需要將四元數單位化
        // We can avoid the expensive sqrt in this case since all 'coefficients' below are products of two q components.
        // That is a square of a square root, so it is possible to avoid that 這是一個平方的平方根,所以可以避免這樣做
        if (length2 != 1.0)
        {
            rlength2 = 2.0/length2;
        }
        else
        {
            rlength2 = 2.0;
        }


        // Source: Gamasutra, Rotating Objects Using Quaternions
        //使用四元數旋轉的演算法
        //http://www.gamasutra.com/features/19980703/quaternions_01.htm


        double wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2;


        // calculate coefficients
        x2 = rlength2*QX;
        y2 = rlength2*QY;
        z2 = rlength2*QZ;


        xx = QX * x2;
        xy = QX * y2;
        xz = QX * z2;


        yy = QY * y2;
        yz = QY * z2;
        zz = QZ * z2;


        wx = QW * x2;
        wy = QW * y2;
        wz = QW * z2;


        // Note.  Gamasutra gets the matrix assignments inverted, resulting
        // in left-handed rotations, which is contrary to OpenGL and OSG's
        // methodology.  The matrix assignment has been altered in the next
        // few lines of code to do the right thing.
        // Don Burns - Oct 13, 2001
        _mat[0][0] = 1.0 - (yy + zz);
        _mat[1][0] = xy - wz;
        _mat[2][0] = xz + wy;




        _mat[0][1] = xy + wz;
        _mat[1][1] = 1.0 - (xx + zz);
        _mat[2][1] = yz - wx;


        _mat[0][2] = xz - wy;
        _mat[1][2] = yz + wx;
        _mat[2][2] = 1.0 - (xx + yy);
    }


#if 0
    _mat[0][3] = 0.0;
    _mat[1][3] = 0.0;
    _mat[2][3] = 0.0;


    _mat[3][0] = 0.0;
    _mat[3][1] = 0.0;
    _mat[3][2] = 0.0;
    _mat[3][3] = 1.0;
#endif
}