3Dプログラミング講座・その3:クォータニオン、四元数



同次座標がw,x,y,z,...の順番で定義されたものの説明です。順番が違うだけで、機能に変わりはありませんが・・・


任意軸の回転としてクォータニオン、四元数を導入する。
これで、記憶容量や演算回数を少なくできるからである。

・基本性質0 クォータニオンQとは

0-1.Q=[w,x,y,z]=w+xi+yj+zk=s+v

と書ける。i,j,kは各軸の虚数単位、sはw、vは[x,y,z]ベクトルである。

虚数単位の性質は

0-2.i^2=j^2=k^2=-1
0-3.ij=-ji=k
0-4.jk=-kj=i
0-5.ki=-ik=j

である。



・基本性質1
・クォータニオンの積

q1=[w1,x1,y1,z1],q2=[w2,x2,y2,z2]のとき

1-1.q1q2=
      (w1w2-x1x2-y1y2-z1z2)
     +(w1x2+x1w2+y1z2-z1y2)i
     +(w1y2-x1z2+y1w2+z1x2)j
     +(w1z2+x1y2-y1x2+z1w2)k

q1=[s1,v1],q2=[v2,v2]のとき

1-2.q1q2=s1s2-v1・v2+s1v2+s2v1+v1×v2



・基本性質2
・共役

2-1.q=[s,v]で与えられたとき、q~=[s,-v]を共役という
2-2.qq~=q~q=q・q=|q|^2=q^2が大きさの2乗である


・基本性質3
・逆元

3-1.逆元q^-1=q~/q^2である



・基本性質4
・回転

4-1.ベクトルP、クォータニオンqにおいて、P.rot(q)=qPq^-1である。

q=[s,v]を単位クォータニオンとすると、q^-1=[s,-v]であり

4-2.
qPq^-1=(s+v)P(s-v)
    =(-v・P+sP+v×P)(s-v)
    =-sv・P+s^2P+s(v×P)+(v・P)v-sPv-(v×P)v
    =s^2P+2s(v×P)+(v・P)v-(v×P×v)

ここで単位ベクトルAとしてv=tAとおけば

4-3.qPq^-1=(s^2-t^2)P+2st(A×P)+2t^2(A・P)A

行列の回転と比較して
t=√((1-cosθ)/2)=sin(θ/2)、またs^2+t^2=1から、s=cos(θ/2)より

4-4.軸Aまわりの角度θのクォータニオンq=cos(θ/2)+Asin(θ/2)

4-5.以上より、まとめると、クォータニオンの変換は

      |1-2y^2-2z^2 2xy-2wz 2xz+2wy|
P.rot(q)=|2yx+2wz 1-2x^2-2z^2 2yz-2wx||P|
      |2zx-2wy 2zy+2wx 1-2x^2-2y^2|

となる。この逆変換は同次行列mとして

4-6.

QSIGN = function(x) {
if(0.0 <= x) return 1.0
return -1.0
}

ret = [( m11 + m22 + m33 + 1.0) / 4.0, ( m.x11 - m22 - m33 + 1.0) / 4.0,
      (-m11 + m22 - m33 + 1.0) / 4.0, (-m11 - m22 + m33 + 1.0) / 4.0)]
if(ret0 < 0) ret0 = 0
if(ret1 < 0) ret1 = 0
if(ret2 < 0) ret2 = 0
if(ret3 < 0) ret3 = 0
ret0 = √(ret0)
ret1 = √(ret1)
ret2 = √(ret2)
ret3 = √(ret3)
r = ret0
if(r < ret1) r = ret1
if(r < ret2) r = ret2
if(r < ret3) r = ret3
if(r == ret0) {
 ret0 *= 1
 ret1 *= QSIGN(m32 - m23)
 ret2 *= QSIGN(m13 - m31)
 ret3 *= QSIGN(m21 - m12)
}
else if(r == ret.q.x[1]) {
 ret0 *= QSIGN(m32 - m23)
 ret1 *= 1
 ret2 *= QSIGN(m21 + m12)
 ret3 *= QSIGN(m13 + m31)
}
else if(r == ret.q.x[2]) {
 ret0 *= QSIGN(m13 - m31)
 ret1 *= QSIGN(m21 + m12)
 ret2 *= 1
 ret3 *= QSIGN(m32 + m23)
}
else {
 ret0 *= QSIGN(m21 - m12)
 ret1 *= QSIGN(m13 + m31)
 ret2 *= QSIGN(m32 + m23)
 ret3 *= 1
}
return ret

である。

4-7.
クォータニオンから回転軸、回転角度を求める

回転角度θは

θ=2 * cos^-1(q.w)

で求められ、回転軸Rは

R=[q.x/sin(θ/2), q.y/sin(θ/2), q.z/sin(θ/2)]

で求められる

θ=0のときはそこから回転軸が求められないので
クォータニオンから回転行列を求め

R=[m32-m23, m13-m31, m21-m12]

これがゼロベクトルであったら

R=[m32, m13, m21]

これがゼロベクトルであったら

R=[m11, m22, m33]

で求められる


4-8.
回転ベクトルで回転演算をするには、クォータニオン演算をそのまま用いることが出来る
回転軸A、θの回転ベクトルのとき、クォータニオンは4-4.から

q=cos(θ/2)+Asin(θ/2)

と表わされ、行列を介さずに、そのクォータニオンで回転演算をすればいいということである


3D回転において
回転ベクトルがあったら
そこからクォータニオンを求めて
回転しちゃえばいいから
行列の出る幕は平行移動だけだよな

平行移動も含めたクォータニオンみたいなのは
出来ないもんかな

多分
θ、3D軸、同次W、3D点の
八元数になると思うけど

これが出来たら4×4同次座標行列は陳腐化するだろうな

思いのほか簡単だった

4-9.
姿勢ベクトルPをこう定義する

同次平行移動T、四元数Qのとき

P=[T Q]

また、この積は

P1=[T1 Q1]、P2=[T2 Q2]、のとき

P1P2=[T1+Q1T2Q1^-1 Q1Q2]

これで、同次座標行列で表現しなくても、それより簡単に計算できるはず


・確かめ算
ポリゴンテスト.htm

ポリゴンテスト.zip

ポリゴンテストFPS・PosVec版.htm

ポリゴンテストFPS・PosVec版.zip


・./javascripts/nas6/testpoly.js


// ……… 前 略 ………

var bp = new N6LVector([1, 0, 0, 8, 1, 0, 0, 0], true);  //pos Box//Box座標

//main loop//メインループ
function GLoop(id){

// ……… 中 略 ………

//姿勢ベクトル演算
  var VecWK = new N6LVector([1,0,0,0,1,0,0,0],true);
  var v = new N6LVector(4, false);
  //unit vector//単位ベクトル
  var ax = new N6LVector(4, true).UnitVec(1);
  var ay = new N6LVector(4, true).UnitVec(2);
  var az = new N6LVector(4, true).UnitVec(3);

  //rot mov obj//物体回転移動
  var q = new N6LQuaternion().UnitQuat().RotAxisQuat(ay, 1.0 * Math.PI / 180.0);
  var pv = new N6LVector([1,0,0,0,q.q.x[0],q.q.x[1],q.q.x[2],q.q.x[3]], true)
  var pvb = new N6LVector([bp.x[0],bp.x[1],bp.x[2],bp.x[3],1,0,0,0], true);
  VecWK = pv.PosVecMul(pvb);  //around y axis rotate 1 degree//y軸回りに1度回転する回転行列を乗算

  //unit vector//単位ベクトル
  ax = ax.UnitVec(1);
  ay = ay.UnitVec(2);
  az = az.UnitVec(3);
  //rot obj//物体回転
  q = new N6LQuaternion().UnitQuat().RotAxisQuat(az, 3.0 * Math.PI / 180.0);
  pv = new N6LVector([1,0,0,0,q.q.x[0],q.q.x[1],q.q.x[2],q.q.x[3]], true)
  pvb = new N6LVector([VecWK.x[0],VecWK.x[1],VecWK.x[2],VecWK.x[3],bp.x[4],bp.x[5],bp.x[6],bp.x[7]], true);
  VecWK = pvb.PosVecMul(pv);  //around z axis rotate 3 degree//z軸回りに3度回転する回転行列を乗算
  q = new N6LQuaternion().UnitQuat().RotAxisQuat(ay, 2.0 * Math.PI / 180.0);
  pv = new N6LVector([1,0,0,0,q.q.x[0],q.q.x[1],q.q.x[2],q.q.x[3]], true)
  VecWK = VecWK.PosVecMul(pv);  //around z axis rotate 3 degree//z軸回りに3度回転する回転行列を乗算
  q = new N6LQuaternion().UnitQuat().RotAxisQuat(ax, 1.0 * Math.PI / 180.0);
  pv = new N6LVector([1,0,0,0,q.q.x[0],q.q.x[1],q.q.x[2],q.q.x[3]], true)
  bp = VecWK.PosVecMul(pv);  //around z axis rotate 3 degree//z軸回りに3度回転する回転行列を乗算


  var tbm = bp.PosVecMatrix();
  var tbx = tbm.Pos();
  var v = tbm.Vector();  //rot vector//回転行列から回転ベクトルを取得

  var angle = tbm.EulerAngle(3, 2, 1);  //rotate order ZYX //回転順番 ZYX

// ↑ 前 半 : ↓ 後 半

  //apply x3dom
  var pos = tbx.ToX3DOM(true);
  var elm = document.getElementById('box0');
  elm.setAttribute('translation', pos.toString());
  var rot = v.ToX3DOM();
  elm = document.getElementById('box1');
  elm.setAttribute('rotation', rot.toString());

  //debug//デバッグ用
  elm = document.getElementById('debug');
  elm.innerText = 
  'EulerAngle(rotate per degree z(3)_y(2)_x(1))\n' + angle.x[0] + ' ' + Math.floor(angle.x[1] * 180.0 / Math.PI) + ' ' + Math.floor(angle.x[2] * 180.0 / Math.PI) + ' ' + Math.floor(angle.x[3] * 180.0 / Math.PI); 

  TMan.timer[id].setalerm(function() { GLoop(id); }, 50);  //reset main loop//メインループ再セット
}

// ……… 後 略 ………


以上で、行列と比べても、正常に動作したから、合っているはず

速さは、項の積の数だけ数えて
行列積が16^2=256で

姿勢ベクトル積が
P1P2=[T1+Q1T2Q1^-1 Q1Q2]

Q1T2Q1^-1
で(四元数→行列)・ベクトル
にT1のベクトル加算と
Q1Q2の四元数の積って計算は

18×3+4^2=70になるけど

本当にそうなって速いか分からん

容量は行列が16個で姿勢ベクトルが8個(同次情報をシェイプすれば7個)


・まとめ

三次元は
3元平行移動T+四元数Qの
7つのパラメタがあれば
完全に記述できる

これを姿勢ベクトルPと定義する
P=[1 T Q]
P1とP2の積、つまり回転は
P1P2=[T1+Q1T2Q1^-1 Q1Q2]
と表記される

速さは、項の積の数だけ数えて
行列積が16^2=256で

姿勢ベクトル積が
P1P2=[T1+Q1T2Q1^-1 Q1Q2]

Q1T2Q1^-1
で(四元数→行列)・ベクトル
にT1のベクトル加算と
Q1Q2の四元数の積って計算は

18×3+4^2=70になる


3D回転テスト・姿勢ベクトルpostest.htm

3D回転テスト・姿勢ベクトルpostest.zip



・基本性質5
・線形補間

5-1.0≦t≦1、初期q1終了q2として

q.lerp(t)=(1-t)q1+tq2



・基本性質6
・球面線形補間

6-1.s=sinθ=√(1-(q1・q2)^2)、0≦t≦1、初期q1終了q2として

q.slerp(t)=(s(1-t)/s)q1+(st/s)q2







3Dプログラミング講座
NAS6LIB
ポリゴンテスト解説・x3dom使い方
その1:ベクトル
その2:行列
その3:クォータニオン
その4:基本総復習・実践
その5:応用その1・メッシュアニメーション、動的テクスチャ
その6:応用その2・GLSL、カスタムシェーダー、キーボード、ファイル
その7:応用その3・ゲームプログラミング、タグの動的追加
その8:応用その4・GLSL、シェーダー、その2

<<prev 行列 : 基本総復習・実践 next>>





戻る