OpenCASCADE Quaternion
eryar@163.com
Abstract. The quaternions are members of a noncommutative division algebra first invented by William Rowan Hamilton. The idea for quaternions occurred to him while he was walking along the Royal Cannal on his way to a meeting of the Irish Academy, and Hamilton was so pleased with his discovery that he scratched the fundamental formula of quaternion algebra. There are several different ways we can express orientation and angular displacement in 3D. Here we discuss the three most important methods-matrices, Euler angles, and quaternions.
Key Words. OpenCASCADE, Quaternion, Euler angles, Rotation, Transformation
1. Introduction
物體在三維空間中經常需要進行一些變換操作，如移動、旋轉等。在CAD軟件中如果提供的便利的交互方式來對模型的位置進行編輯，則軟件的用戶體驗就很好。對模型位置編輯的結果需要在計算機中保存起來，這就需要用一種方式來記錄模型的變換。如下圖所示：
Figure 1.1 Modify the location of a Valve
上圖所示為三維工廠軟件PDMS中使用交互的方式來修改一個閥件的位置，直接使用鼠標拖動高亮的箭頭即可實現對閥件位置的編輯。
物體在三維空間中的旋轉變換操作常見的有三種表示方法：Matrix、Euler Angles、Quaternion。每種表示方式各用利弊，根據需要選擇合適的旋轉的表示方式。本文詳細介紹了這三種方式及在OpenCASCADE中使用Quaternion和在這三種表示方式之間進行轉換。
2.Rotation in Matrix Form
在3D空間中描述坐標系的方位（orientation）的一種方法就是列出這個坐標系的基向量，這些基向量是用其他的坐標系來描述的。用這些基向量構成一個3X3矩陣，然后就能用矩陣的形式來描述方位了。換句話說，能用一個旋轉矩陣來描述這兩個坐標系之間的相對方位，這個旋轉矩陣用于把一個坐標系中的向量變換到另一個坐標系中，如下圖所示：
Figure 2.1 Defining an orientation using a matrix
Figure 2.2 A Valve Orientation in PDMS
由上圖可知，在PDMS中對模型的方位的保存也是采用了矩陣形式，其中X is E，Y is N 60U， Z is S 30 U這其實是三個向量。下面給出繞任意軸旋轉一定角度的矩陣表示的證明：
Figure 2.3 Rotating a vector about an arbitrary axis
根據向量的運算規則容易推出：
則繞軸旋轉的矩陣形式如下：
其中在OpenCASCADE的類gp_Mat中實現代碼如下所示：
void gp_Mat::SetRotation (const gp_XYZ& Axis,
const Standard_Real Ang)
{
// Rot = I + sin(Ang) * M + (1. - cos(Ang)) * M*M
// avec M . XYZ = Axis ^ XYZ
// const Standard_Address M = (Standard_Address)&(matrix[0][0]);
gp_XYZ V = Axis.Normalized();
SetCross (V);
Multiply (sin(Ang));
gp_Mat Temp;
Temp.SetScale (1.0);
Add (Temp);
Standard_Real A = V.X();
Standard_Real B = V.Y();
Standard_Real C = V.Z();
Temp.SetRow (1, gp_XYZ(- C*C - B*B, A*B, A*C ));
Temp.SetRow (2, gp_XYZ( A*B, -A*A - C*C, B*C ));
Temp.SetRow (3, gp_XYZ( A*C, B*C, - A*A - B*B));
Temp.Multiply (1.0 - cos(Ang));
Add (Temp);
}
3.Rotation with Euler Angles
用Euler角的方式來表示旋轉這項技術是以著名數學家Leonhard Euler(1707~1783)來命名的，他證明了角位移序列等價于單個角位移，即可以用一個合成的變換來表示多個連續的變換，證明過程詳見蘇步青《應用幾何教程》。
Euler Angles的基本思想是將角位移分解為繞三個互相垂直的三個旋轉組成的序列。這聽起來有點復雜，其實是非常直觀的，這也正是Euler Angle易于使用的優點之一。Euler Angle將方位Orientation分解為繞三個垂直軸的旋轉，那么是哪三個軸？按什么順序？其實任意三個軸和任意順序都可以，但最有意義的是使用笛卡爾坐標系按一定順序所組成的旋轉序列。最常用的約定是所謂的“heading-picth-bank”，如下圖所示為給定heading，picth和bank角度后，可以用四步法來確定Euler Angle對應的Orientation：
Figure 3.1 Step 1: An object in its identity orientation
Figure 3.2 Step 2: Heading is the first rotation and rotates about the vertical axis(y-axis)
Figure 3.3 Step 3: Pitch is the second rotation and rotates about the object laterial axis(x-axis)
Figure 3.4 Step 4: Bank is the third and rotates about the object longitudinal axis(z-axis)
heading-pitch-bank系統不是唯一的Euler Angle系統。繞任意三個互相垂直的任意旋轉序列都能定義一個方位orientation。所以多種選擇導致了Euler Angle約定的多樣性。如常用的術語roll-pitch-yaw，其中roll等價于bank，yaw基本上等價于heading，他的順序與heading-pitch-bank相反。
因為Euler Angle的易用性，只需要約定旋轉序列和三個角度即可表示方位了。可以仿照上述變換過程應用Euler Angle來實現模型旋轉編輯的交互操作，實現交互方式友好的操作。即當鼠標移動到高亮的旋轉handle上時，就可以繞一個軸旋轉一定角度，不一定是按heading-pitch-bank的順序來旋轉。如下圖所示：
Figure 5. Model Editor in PDMS
Figure 6. Euler Angle in Model Editor
由上圖可知，對模型進行旋轉時，實時顯示的角度正是Euler Angle。如果約定了Euler Angle的順序，如為heading-pitch-bank，只需要三個實數即可表示方位orientation，這在將方位數據保存到文件時有很大優勢，可以節省大量存儲空間。如下圖所示：
Figure 7. Orientation Properties in AVEVA Plant/PDMS
由上圖可知，PDMS的數據庫中保存orientation的方式使用了Euler Angle的方式，這樣與矩陣的方式相比，數據量要少兩倍，因為一般的模型都有這個屬性，這樣下來，節省的存儲空間比較可觀。
因為Euler Angles的方式最簡單直觀，就是三個角度，所以一般都是將Euler Angle轉換到其他的形式，如Matrix和Quaternion。
4.Quaternions
一個Quaternion包含一個標量分量和一個3D向量分量，一般記標量分量為w，向量分量為V或分開的x,y,z，如Q=[w, V]或Q=[w,(x,y,z)]。
Quaternion能被解釋為角位移的軸－角對方式。然而，旋轉軸和角度不是直接存儲在Quaternion的四個數中，它們的確在Quaternion中，但是不是那么直接，其關系式為：
根據這個關系式即可實現繞軸旋轉一定角度的Quaternion，其中OpenCASCADE中的實現代碼如下所示：
//=======================================================================
//function : SetVectorAndAngle
//purpose :
//=======================================================================
void gp_Quaternion::SetVectorAndAngle (const gp_Vec& theAxis,
const Standard_Real theAngle)
{
gp_Vec anAxis = theAxis.Normalized();
Standard_Real anAngleHalf = 0.5 * theAngle;
Standard_Real sin_a = Sin (anAngleHalf);
Set (anAxis.X() * sin_a, anAxis.Y() * sin_a, anAxis.Z() * sin_a, Cos (anAngleHalf));
}
為了將Quaternion轉換到矩陣形式，可以利用繞任意軸旋轉的矩陣，公式如下所示：
矩陣是由旋轉軸n和旋轉角度θ表示的，但是Quaternion是由下述分量表示的：
將R(n,θ)變換到R(w,x,y,z)是一個技巧性很強的推導，如果只是為了使用矩陣，那么就不必理解矩陣是如何推導的。如果對推導過程感興趣，可參考Fletcher Dunn,Ian Parberry.所著《3D Math Primer for Graphics and Game Development》，其中有詳細推導說明。下面就直接給出推導結果：
Figure 4.1 Converting a quaternion to a 3x3 matrix
在OpenCASCADE中取得Quaternion的矩陣數據是通過函數GetMatrix，其實現代碼如下所示：
//=======================================================================
//function : GetMatrix
//purpose :
//=======================================================================
gp_Mat gp_Quaternion::GetMatrix () const
{
Standard_Real wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2;
Standard_Real s = 2.0 / SquareNorm();
x2 = x * s; y2 = y * s; z2 = z * s;
xx = x * x2; xy = x * y2; xz = x * z2;
yy = y * y2; yz = y * z2; zz = z * z2;
wx = w * x2; wy = w * y2; wz = w * z2;
gp_Mat aMat;
aMat (1, 1) = 1.0 - (yy + zz);
aMat (1, 2) = xy - wz;
aMat (1, 3) = xz + wy;
aMat (2, 1) = xy + wz;
aMat (2, 2) = 1.0 - (xx + zz);
aMat (2, 3) = yz - wx;
aMat (3, 1) = xz - wy;
aMat (3, 2) = yz + wx;
aMat (3, 3) = 1.0 - (xx + yy);
// 1 division 16 multiplications 15 addidtions 12 variables
return aMat;
}
由代碼可知，即是對上述轉換公式的直接實現。如其注釋中所言，轉換涉及到了1個除法、16個乘法、15個加法和12個變量。
OpenCASCADE中還實現了從Euler Angles轉換到Quaternion的功能，實現代碼如下所示：
//=======================================================================
//function : GetEulerAngles
//purpose :
//=======================================================================
void gp_Quaternion::GetEulerAngles (const gp_EulerSequence theOrder,
Standard_Real& theAlpha,
Standard_Real& theBeta,
Standard_Real& theGamma) const
{
gp_Mat M = GetMatrix();
gp_EulerSequence_Parameters o = translateEulerSequence (theOrder);
if ( o.isTwoAxes )
{
double sy = sqrt (M(o.i, o.j) * M(o.i, o.j) + M(o.i, o.k) * M(o.i, o.k));
if (sy > 16 * DBL_EPSILON)
{
theAlpha = ATan2 (M(o.i, o.j), M(o.i, o.k));
theGamma = ATan2 (M(o.j, o.i), -M(o.k, o.i));
}
else
{
theAlpha = ATan2 (-M(o.j, o.k), M(o.j, o.j));
theGamma = 0.;
}
theBeta = ATan2 (sy, M(o.i, o.i));
}
else
{
double cy = sqrt (M(o.i, o.i) * M(o.i, o.i) + M(o.j, o.i) * M(o.j, o.i));
if (cy > 16 * DBL_EPSILON)
{
theAlpha = ATan2 (M(o.k, o.j), M(o.k, o.k));
theGamma = ATan2 (M(o.j, o.i), M(o.i, o.i));
}
else
{
theAlpha = ATan2 (-M(o.j, o.k), M(o.j, o.j));
theGamma = 0.;
}
theBeta = ATan2 (-M(o.k, o.i), cy);
}
if ( o.isOdd )
{
theAlpha = -theAlpha;
theBeta = -theBeta;
theGamma = -theGamma;
}
if ( ! o.isExtrinsic )
{
Standard_Real aFirst = theAlpha;
theAlpha = theGamma;
theGamma = aFirst;
}
}
下面給出在OpenCASCADE中使用Quaternion的具體示例，代碼如下所示：
/*
* Copyright (c) 2013 to current year. All Rights Reserved.
*
* File : Main.cpp
* Author : eryar@163.com
* Date : 2014-11-29 10:18
* Version : OpenCASCADE6.8.0
*
* Description : Test OpenCASCADE quaternion.
*
* Key Words : OpenCASCADE, Quaternion
*
*/
#define WNT
#include <gp_Quaternion.hxx>
#pragma comment(lib, "TKernel.lib")
#pragma comment(lib, "TKMath.lib")
void TestQuaternion(void)
{
gp_Quaternion aQuaternion;
// create quaternion by axis-angle.
aQuaternion.SetVectorAndAngle(gp_Vec(1.0, 0.0, 0.0), M_PI_2);
// convert quaternion to matrix.
gp_Mat aMatrix = aQuaternion.GetMatrix();
Standard_Real aYaw = 0.0;
Standard_Real aPitch = 0.0;
Standard_Real aRoll = 0.0;
// convert quaternion to Euler Angles.
aQuaternion.GetEulerAngles(gp_YawPitchRoll, aYaw, aPitch, aRoll);
}
int main(int argc, char* argv[])
{
TestQuaternion();
return 0;
}
5.Conclusions
綜上所述，Euler Angles最容易被使用，當需要為世界中的物體指定方位時，Euler Angles能大簡化人機交互，包括直接的鍵盤、鼠標輸入及在調試中測試。
如果需要要坐標系之間進行轉換向量，那么就選矩陣形式。當然這并不意味著你不能使用其他格式來保存方位，并在需要的時候轉換到矩陣形式。另一種方法是用Euler Angles作為方位的“主拷貝”，并同時維護一個旋轉矩陣，當Euler Angles發生變化時矩陣也同時進行更新。
當需要大量保存方位數據時，就使用Euler Angles或Quaternion。Euler Angles將少占用25%的空間，但它在轉換到矩陣時要稍微慢點。如果動畫數據需要進行坐標系之間的連接，Quaternion可能是最好的選擇了。
平滑的插值只能用Quaternion來完成，這在OpenSceneGraph中有大量的應用。如果你使用其他形式，也可先轉換成Quaternion再進行插值，插值完成后再轉換回原來的形式。
OpenCASCADE的Quaternion類中實現了Matrix, Euler Angles的轉換，即通過gp_Quaternion即可將旋轉操作在這三種形式之間進行轉換。
6. References
1. WolframMathWorld, http://mathworld.wolfram.com/Quaternion.html
2. Ken Shoemake. Conversion between quaternion and Euler angles. Graphics Gems IV, P222-22
http://tog.acm.org/resources/GraphicsGems/gemsiv/euler_angle/EulerAngles.c
3. 蘇步青, 華宣積. 應用幾何教程. 復旦大學出版計. 2012
4. 丘維聲. 解析幾何. 北京大學出版社. 1996
5. 同濟大學應用數學系編. 線性代數（第四版）. 高等教育出版社. 2003
6. Fletcher Dunn,Ian Parberry. 3D Math Primer for Graphics and Game Development. CRC Press
7. 史銀雪,陳洪,王榮靜. 3D數學基礎：圖形與游戲開發. 清華大學出版社. 2005
8. 蘇步青. 神奇的符號. 湖南少年兒童出版社. 2010