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


クォータニオン、四元数とは本質的には何か?



軸の直交と虚軸関係

まずこれを見てもらって

反転軸がX^2=-I
直交軸がX=iI

という理解で四元数を考えると
wxyzの各パラメタがお互い直交関係であるということです
したがってユークリッド系と親和性が高く
ユークリッド系の回転変換に使えます







 ■■■ 3Dプログラミング入門講座 ■■■ 
NAS6LIB
ポリゴンテスト解説・x3dom使い方
その1:ベクトル
その2:行列
その3:クォータニオン
その4:基本総復習・実践
その5:応用その1・メッシュアニメーション、動的テクスチャ
その6:応用その2・GLSL、カスタムシェーダー、キーボード、ファイル
その7:応用その3・ゲームプログラミング、タグの動的追加
その8:応用その4・GLSL、シェーダー、その2
その9:物理演算その1・電卓で相対性理論を解く方法
その10:物理演算その2・相対性理論的ニュートン力学
その11:物理演算その3・ケプラー方程式で惑星軌道シミュレーターを作る

その12:物理演算その4・ルンゲクッタ法で作った 相対性理論的ニュートン力学物理エンジンで惑星軌道シミュレーターを作る

その13:経路探索(A*:A-STAR)&巡回セールスマン問題 : 巨大サイズ : くろにゃんこ

その14:プログラミングにおける配列テーブルテクニック
その15:javascriptのクラス活用法
その16:透視射影公式テスト

その17:ケプラー方程式カプセルライブラリ使用法
その18:CSVファイル処理
その19:物理演算その5・重力多体問題
その20:同次座標について(3D座標系の基本の基本)
その21:おさらいコモンクラスの宣言
その22:物理エンジンライブラリ解説(ケプラー方程式・ルンゲクッタ・相対論的万有引力)


 ■■■ THREE.JSプログラミング講座 ■■■ 
THREE.JSテスト解説・THREE.JS使い方
THREE.JS examplesをいじってみた(フレネル反射透過シェーダー)

THREE.JS (半透明シェーダー)

THREE.JS 3D演算で必要な計算(具体例)★とても重要★
THREE.JS THREE-VRM をいじってみた



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





戻る