分類： VC++備忘錄 Graphics 2013-05-20 20:25 14396人閱讀 評論(1) 收藏 舉報
Quaternion(四元數)
Quaternion 的定義
四元數一般定義如下：
q=w+xi+yj+zk
其中 w,x,y,z是實數。同時，有:
i*i=-1
j*j=-1
k*k=-1
四元數也可以表示為：
q=[w,v]
其中v=(x,y,z)是矢量，w是标量，雖然v是矢量，但不能簡單的了解為3D空間的矢量，它是4維空間中的的矢量，也是非常不容易想像的。
通俗的講，一個四元數（Quaternion）描述了一個旋轉軸和一個旋轉角度。這個旋轉軸和這個角度可以通過 Quaternion::ToAngleAxis轉換得到。當然也可以随意指定一個角度一個旋轉軸來構造一個Quaternion。這個角度是相對于機關四元數而言的，也可以說是相對于物體的初始方向而言的。
當用一個四元數乘以一個向量時，實際上就是讓該向量圍繞着這個四元數所描述的旋轉軸，轉動這個四元數所描述的角度而得到的向量。
四元組的優點
有多種方式可表示旋轉，如 axis/angle、歐拉角(Euler angles)、矩陣(matrix)、四元組等。 相對于其它方法，四元組有其本身的優點：
四元數不會有歐拉角存在的 gimbal lock 問題
四元數由4個數組成，旋轉矩陣需要9個數
兩個四元數之間更容易插值
四元數、矩陣在多次運算後會積攢誤差，需要分别對其做規範化(normalize)和正交化(orthogonalize)，對四元數規範化更容易
與旋轉矩陣類似，兩個四元組相乘可表示兩次旋轉
Quaternion 的基本運算
Normalizing a quaternion
// normalising a quaternion works similar to a vector. This method will not do anything
// if the quaternion is close enough to being unit-length. define TOLERANCE as something
// small like 0.00001f to get accurate results
void Quaternion::normalise()
{
// Don't normalize if we don't have to
float mag2 = w * w + x * x + y * y + z * z;
if ( mag2!=0.f && (fabs(mag2 - 1.0f) > TOLERANCE)) {
float mag = sqrt(mag2);
w /= mag;
x /= mag;
y /= mag;
z /= mag;
}
}
The complex conjugate of a quaternion
// We need to get the inverse of a quaternion to properly apply a quaternion-rotation to a vector
// The conjugate of a quaternion is the same as the inverse, as long as the quaternion is unit-length
Quaternion Quaternion::getConjugate()
{
return Quaternion(-x, -y, -z, w);
}
Multiplying quaternions
// Multiplying q1 with q2 applies the rotation q2 to q1
Quaternion Quaternion::operator* (const Quaternion &rq) const
{
// the constructor takes its arguments as (x, y, z, w)
return Quaternion(w * rq.x + x * rq.w + y * rq.z - z * rq.y,
w * rq.y + y * rq.w + z * rq.x - x * rq.z,
w * rq.z + z * rq.w + x * rq.y - y * rq.x,
w * rq.w - x * rq.x - y * rq.y - z * rq.z);
}
Rotating vectors
// Multiplying a quaternion q with a vector v applies the q-rotation to v
Vector3 Quaternion::operator* (const Vector3 &vec) const
{
Vector3 vn(vec);
vn.normalise();
Quaternion vecQuat, resQuat;
vecQuat.x = vn.x;
vecQuat.y = vn.y;
vecQuat.z = vn.z;
vecQuat.w = 0.0f;
resQuat = vecQuat * getConjugate();
resQuat = *this * resQuat;
return (Vector3(resQuat.x, resQuat.y, resQuat.z));
}
How to convert to/from quaternions1
Quaternion from axis-angle
// Convert from Axis Angle
void Quaternion::FromAxis(const Vector3 &v, float angle)
{
float sinAngle;
angle *= 0.5f;
Vector3 vn(v);
vn.normalise();
sinAngle = sin(angle);
x = (vn.x * sinAngle);
y = (vn.y * sinAngle);
z = (vn.z * sinAngle);
w = cos(angle);
}
Quaternion from Euler angles
// Convert from Euler Angles
void Quaternion::FromEuler(float pitch, float yaw, float roll)
{
// Basically we create 3 Quaternions, one for pitch, one for yaw, one for roll
// and multiply those together.
// the calculation below does the same, just shorter
float p = pitch * PIOVER180 / 2.0;
float y = yaw * PIOVER180 / 2.0;
float r = roll * PIOVER180 / 2.0;
float sinp = sin(p);
float siny = sin(y);
float sinr = sin(r);
float cosp = cos(p);
float cosy = cos(y);
float cosr = cos(r);
this->x = sinr * cosp * cosy - cosr * sinp * siny;
this->y = cosr * sinp * cosy + sinr * cosp * siny;
this->z = cosr * cosp * siny - sinr * sinp * cosy;
this->w = cosr * cosp * cosy + sinr * sinp * siny;
normalise();
}
Quaternion to Matrix
// Convert to Matrix
Matrix4 Quaternion::getMatrix() const
{
float x2 = x * x;
float y2 = y * y;
float z2 = z * z;
float xy = x * y;
float xz = x * z;
float yz = y * z;
float wx = w * x;
float wy = w * y;
float wz = w * z;
// This calculation would be a lot more complicated for non-unit length quaternions
// Note: The constructor of Matrix4 expects the Matrix in column-major format like expected by
// OpenGL
return Matrix4( 1.0f - 2.0f * (y2 + z2), 2.0f * (xy - wz), 2.0f * (xz + wy), 0.0f,
2.0f * (xy + wz), 1.0f - 2.0f * (x2 + z2), 2.0f * (yz - wx), 0.0f,
2.0f * (xz - wy), 2.0f * (yz + wx), 1.0f - 2.0f * (x2 + y2), 0.0f,
0.0f, 0.0f, 0.0f, 1.0f)
}
Quaternion to axis-angle
// Convert to Axis/Angles
void Quaternion::getAxisAngle(Vector3 *axis, float *angle)
{
float scale = sqrt(x * x + y * y + z * z);
axis->x = x / scale;
axis->y = y / scale;
axis->z = z / scale;
*angle = acos(w) * 2.0f;
}
Quaternion 插值
線性插值
最簡單的插值算法就是線性插值，公式如：
q(t)=(1-t)q1 + t q2
但這個結果是需要規格化的，否則q(t)的機關長度會發生變化，是以
q(t)=(1-t)q1+t q2 / || (1-t)q1+t q2 ||
球形線性插值
盡管線性插值很有效，但不能以恒定的速率描述q1到q2之間的曲線，這也是其弊端，我們需要找到一種插值方法使得q1->q(t)之間的夾角θ是線性的，即θ(t)=(1-t)θ1+t*θ2，這樣我們得到了球形線性插值函數q(t)，如下：
q(t)=q1 * sinθ(1-t)/sinθ + q2 * sinθt/sineθ
如果使用D3D，可以直接使用 D3DXQuaternionSlerp 函數就可以完成這個插值過程。
用 Quaternion 實作 Camera 旋轉
總體來講，Camera 的操作可分為如下幾類：
沿直線移動
圍繞某軸自轉
圍繞某軸公轉
下面是一個使用了 Quaternion 的 Camera 類：
class Camera {
private:
Quaternion m_orientation;
public:
void rotate (const Quaternion& q);
void rotate(const Vector3& axis, const Radian& angle);
void roll (const GLfloat angle);
void yaw (const GLfloat angle);
void pitch (const GLfloat angle);
};
void Camera::rotate(const Quaternion& q)
{
// Note the order of the mult, i.e. q comes after
m_Orientation = q * m_Orientation;
}
void Camera::rotate(const Vector3& axis, const Radian& angle)
{
Quaternion q;
q.FromAngleAxis(angle,axis);
rotate(q);
}
void Camera::roll (const GLfloat angle) //in radian
{
Vector3 zAxis = m_Orientation * Vector3::UNIT_Z;
rotate(zAxis, angleInRadian);
}
void Camera::yaw (const GLfloat angle) //in degree
{
Vector3 yAxis;
{
// Rotate around local Y axis
yAxis = m_Orientation * Vector3::UNIT_Y;
}
rotate(yAxis, angleInRadian);
}
void Camera::pitch (const GLfloat angle) //in radian
{
Vector3 xAxis = m_Orientation * Vector3::UNIT_X;
rotate(xAxis, angleInRadian);
}
void Camera::gluLookAt() {
GLfloat m[4][4];
identf(&m[0][0]);
m_Orientation.createMatrix (&m[0][0]);
glMultMatrixf(&m[0][0]);
glTranslatef(-m_eyex, -m_eyey, -m_eyez);
}
用 Quaternion 實作 trackball
用滑鼠拖動物體在三維空間裡旋轉，一般設計一個 trackball，其内部實作也常用四元數。
class TrackBall
{
public:
TrackBall();
void push(const QPointF& p);
void move(const QPointF& p);
void release(const QPointF& p);
QQuaternion rotation() const;
private:
QQuaternion m_rotation;
QVector3D m_axis;
float m_angularVelocity;
QPointF m_lastPos;
};
void TrackBall::move(const QPointF& p)
{
if (!m_pressed)
return;
QVector3D lastPos3D = QVector3D(m_lastPos.x(), m_lastPos.y(), 0.0f);
float sqrZ = 1 - QVector3D::dotProduct(lastPos3D, lastPos3D);
if (sqrZ > 0)
lastPos3D.setZ(sqrt(sqrZ));
else
lastPos3D.normalize();
QVector3D currentPos3D = QVector3D(p.x(), p.y(), 0.0f);
sqrZ = 1 - QVector3D::dotProduct(currentPos3D, currentPos3D);
if (sqrZ > 0)
currentPos3D.setZ(sqrt(sqrZ));
else
currentPos3D.normalize();
m_axis = QVector3D::crossProduct(lastPos3D, currentPos3D);
float angle = 180 / PI * asin(sqrt(QVector3D::dotProduct(m_axis, m_axis)));
m_axis.normalize();
m_rotation = QQuaternion::fromAxisAndAngle(m_axis, angle) * m_rotation;
m_lastPos = p;
}
---------------------------------------------------------------------------------------------------------
每一個機關四元數都可以對應到一個旋轉矩陣
機關四元數q=（s，V）的共轭為q*=（s，-V）
機關四元數的模為||q||=1；
四元數q=(s,V)的逆q^(-1)=q*/(||q||)=q*
一個向量r，沿着向量n旋轉a角度之後的向量是哪個（假設為v），這個用四元數可以輕松搞定
構造兩個四元數q=（cos(a/2),sin(a/2)*n）,p=(0,r)
p`=q * p * q^(-1) 這個可以保證求出來的p`也是（0，r`）形式的，求出的r`就是r旋轉後的向量
另外其實對p做q * p * q^(-1)操作就是相當于對p乘了一個旋轉矩陣，這裡先假設 q=（cos(a/2),sin(a/2)*n）=（s，（x, y, z））
兩個四元數相乘也表示一個旋轉
Q1 * Q2 表示先以Q2旋轉，再以Q1旋轉
則這個矩陣為
同理一個旋轉矩陣也可以轉換為一個四元數，即給你一個旋轉矩陣可以求出（s，x，y，z）這個四元數，
方法是：
給定任意機關軸q(q1,q2,q3)(向量),求向量p(x,y,z)(或點p)饒q旋轉theta角度的變換後的新向量p'(或點p')：
1.用四元數工具：
-------------------------------------------------------------------------
結論：構造四元數變換p'= q*p*q-1，（p,q是由向量p,q擴充成的四元數）。那麼，p'轉換至對應的向量(或點)就是變換後的新向量p'(或點p')。
其中，p',q,p,q-1均為四元數。q由向量q擴充,為q=(cos(theta/2),sin(theta/2)*q)，p由向量p擴充,為p=(0,x,y,z),q-1為q的逆，因為q為機關四元數，是以q-1=q*=(cos(theta/2),-sin(theta/2)*q)。
http://www.linuxgraphics.cn/opengl/opengl_quaternion.html
http://blog.csdn.net/qq960885333/article/details/8191448
http://blog.csdn.net/jiexuan357/article/details/7727634