せっかくだから剛体回転でも、と思って準備

MSDNのクオータニオン関係の説明が詳細を欠いていて、係数の役割や回転変換との関係付けなど細かな点で不明な所があるので、剛体の回転も含めて整理しようとしたのだが、やっぱり集中力が阻害されて本をまともに読めない。まあ復習だしと思って自分ではじめからやってみる。

まず、三次元回転行列とは三次元の基底ベクトルを並べて作る。つまり、その行または列は正規直交する。従って回転行列RはRR^T=Iで特徴付けられる。

単位行列Iの近くの回転行列だけを考える。つまり、式(I+εM)(I+εM)^T=Iが、εを小さくしていった時にどんな条件式と等価になるかを考える。それは、
M=-M^T
である。この条件式はMの定数倍に対して保存される。つまり、Mに対して成り立てばεMに対しても成り立つ。つまり、単位行列Iの``近く''の回転行列の集合は
{I+M, where M=-M^T, Mは``小さい''}
と書ける。そこで、線形空間
L={M, where M=-M^T}
を考える事が重要になる。

反対称という事は、対角成分は0で上三角部分を決めると下三角部分が決まるから自由に決められるパラメタは三つしかない。これは、この線形空間が三次元空間である事を意味する。なんて勿体つけなくても、基底は直ぐに見つかって、

[0, 0,0]
[0, 0,1]
[0,-1,0]

[ 0,0,1]
[ 0,0,0]
[-1,0,0]

[ 0,1,0]
[-1,0,0]
[ 0,0,0]

がLを張る。

x,y,zが対称的に扱われるように、以下、右手系における右ネジ回転を基準とする事にする。これでこのエントリにおける不明瞭な点はなくなるだろう。後でWPFのクオータニオンと対応付ける必要はあるが。また、点の座標の変換と座標系の変換を明確に区別する。

(x,y,z)にある点をx軸でθ回転した時移る先の座標(x',y',z')は

[x'] [1,0  　,0   　]
[y']=[0,cosθ,-sinθ]
[z'] [0,sinθ,cosθ ]

と表される。この事から、x軸中心の微小なθ回転は

  [0,0,0]
I+[0,0,-θ]
  [0,θ,0]

と表される。他の軸中心の微小回転も同様に求められる。そこでLの基底を次のように取り直すと便利だ。

   [0,0, 0]
Sx=[0,0,-1]
   [0,1, 0]

   [ 0,0,1]
Sy=[ 0,0,0]
   [-1,0,0]

   [0,-1,0]
Sz=[1, 0,0]
   [0, 0,0]

この場合、積はIを加えれば閉じていて、反対称である。つまり
SxSx=-I, SxSy=Sz, SySx=-Sz, ... .

恒等変換の``近く''の回転変換を次々と作用させていく事で、恒等変換から``遠く''の回転変換も作る事ができる。剛体の回転運動はこのようにイメージできる。対応する概念は単位行列の``近く''の回転行列の複数回の積だ。つまり、x(t), y(t), z(t)を与えて、
Π(I+δt(x(t_n)Sx+y(t_n)Sy+z(t_n)S_z)), where δt=t_n-t_(n-1) ...(1)
によって単位行列から``連続に''別の回転行列に移る事ができる。ここでSx, S_y, S_zが単位ラジアンあたりの各軸周りの回転に対応していた事を思い出せば、x(t), y(t), z(t)は角速度ベクトルの各成分を与えると分かる。

この計算を実行するのは難しいが、ひとつだけ分かる事は、I, Sx, Sy, Szの積が閉じているために、この結果は必ず、I, Sx, Sy, Szの線形結合になるという事だ(Lの基底と違い、Iが加わる事に注意！)。これは重要な事だ。というのも、回転行列の積を繰り返すと計算誤差により徐々に回転行列の条件＝正規直交性を満たさなくなってしまうからだ。しかし、これら四つの基底に関する係数成分を時々刻々と計算するならば、計算誤差がどれだけあろうとも、それは回転変換に厳密な意味で対応し続ける！こうして四元数が登場する。

少なくとも剛体のシミュレーションに関する限り式(1)の計算結果を具体的に表す方法がない事はほとんど問題にならない。なぜならば、実際に、その時点での回転行列に恒等変換の近くの回転行列を乗じて、回転変換を更新する処理を繰り返せば良いからだ。変数に更新量を加算したり乗じたりしてループを回す類の処理は、シミュレーションにおける最も基本的なアルゴリズムで誰もが良く慣れているはずだ。また、すでに注意したように、回転行列自身の代わりに基底行列I, Sx, Sy, Szの係数を保存し、回転行列の積を計算する時には、基底行列の間の積の関係を使って係数の組から結果の係数を求めるようにすれば、係数に対応する回転行列は厳密に何らかの回転変換に対応し続ける。

存在する問題は二つだ。一つは、(1)式の積に現れる各項は近似的にしか回転行列に対応しないという事。もう一つは初期姿勢を決めるために(1)式の極限の具体的な結果が一つ欲しいという事。二つ目が解決されれば、その結果は、一つ目の問題の解決にも貢献するだろう。つまり、本物の微小な回転変換を繰り返し合成して一般の回転変換を作る事ができるようになる。

さて、(1)式の計算結果を具体的に求められる簡単なケースが一つある。それは、x(t), y(t), z(t)が一定である場合で、結果は行列の指数関数になる。その意味も明確だ。これは、角速度ベクトル[x, y, z]による、軸・角速度一定の回転を表す。特に、|[x,y,z]|=1ならば、時刻t=θにおいて角度θだけ回転した事になる。これがいわゆるクオータニオンと、軸／回転角で表した回転の関係だ。

ふう。次回は、せっかくだからSU(2)でやらずにこのまま気合で対角化して指数関数を求めるか。

ちょ、

足が、足がやばい・・・

クオータニオンを使って姿勢をランダムにした

クオータニオンで一様な回転を生成するには、単にその四つの係数を四次元球面上に一様に分布させれば良い(⇒これは誤り。回転軸を球面上一様分布で、回転角を一様分布で生成してクオータニオンを生成すべき。だったら、AxisAngleRoatate3Dで良かったんだな。まあ良いや。回転軸は確かに一様だし)。実際には、四つの係数を[0,1]の一様分布で生成する事を長さが1以下になるまで繰り返して、その後まずクオータニオンを作ってしまってからNormalize()を呼べば良い。

というわけで今現在のコードで、ミクさんのモデルデータは重心が原点にくるように移動させてあって、それをContentsとするVisual3DのTransformにTransform3DGroupを設定して、そのChildrenをRotateTransform3D⇒TranslateTransform3Dの順に設定し、この二つの変換を支配するQuaternionおよびPoint3Dを公開して、UpdatePose()を呼ぶと変換に反映するようにした。現在はPoint3Dのみを重力加速度に従って動かしているけれど、このデータ構造なら、オイラーの剛体方程式に従って剛体回転も容易に扱えるはず。やる気はしないけれど。

コードと実行ファイルはARWithWPFMikuEmitted.zipを詳細サイトのダウンロードページから取ってきてください。

ちなみに、この画像ではボタンの``Emit Miku!''の綴りが``Emmit Miku!''となっていて格好悪いけれど、アップしたソースでは修正してあります。

ARToolkitPlusで体感ゲームを作る予定は未定だけど準備だけ

ミクさんがマーカの位置からランダムな初速度を与えられて射出されて重力で自由落下するぞ。

ミクさんは三次元CG@七葉の4169さんのローポリなデータをお借りしています。

しかしWPFのStoryboard使わずに普通にオンフレームアニメーションにしてしまったので、プログラム的にはあまり面白くない。一連のARToolkit+WPFの実験はWPFを活かす事に主眼があるので続くかどうかは未定。

演算子とジェネリック

演算子の実装を表す型制約がないから何かやだねー、というお話。2chのC#板で出ていた。
この手の話は、C++風にダックタイピングが出来たら・・・という話にもなるのだが、そういう方針で、まあマクロのお化け的なC++のテンプレートにおけるダックタイピングとは違くてつまりは実行時に動的にメソッドを生成するという方針だけど、どなたかがExpression TreeのCompileメソッドを使うやり方を示していた。

で、それならCodeDOMでもできるじゃん、と思って書いたのがこんなコード。

    class Multiplier<T>
    {
        public T Mult(T l, T r)
        {
            string nameOfT = typeof(T).FullName;
            string assemblyOfT = typeof(T).Assembly.ManifestModule.FullyQualifiedName;
            CompilerResults cr = null; 

            using (CodeDomProvider provider
                = new CSharpCodeProvider(new Dictionary() { { "CompilerVersion", "v3.5" } })) {
                CompilerParameters cp = new CompilerParameters() {
                    GenerateInMemory = true
                };
                cp.ReferencedAssemblies.AddRange(new string[] { "System.Core.dll", "mscorlib.dll", assemblyOfT });
                string className = "Multiplier_" + nameOfT.Replace('.', '_');
                string namespaceName="GenericTest";
                string qualifiedClassName = namespaceName + "." + className;
                string source = "namespace " +namespaceName+"{ public class " + className + "{public " + nameOfT + " Mult(" + nameOfT + " l, " + nameOfT + @" r){return l*r;}}}";

                cr = provider.CompileAssemblyFromSource(cp, source);
                if (cr.Errors.HasErrors) {
                    throw new Exception(String.Join(Environment.NewLine, cr.Errors.OfType().Select(e => e.ToString()).ToArray()));
                }
                Type typeOfMultiplier = cr.CompiledAssembly.GetType(qualifiedClassName);
                Object multiplier = Activator.CreateInstance(typeOfMultiplier);
                MethodInfo multiplierMethodMult = typeOfMultiplier.GetMethod("Mult");
                return (T)multiplierMethodMult.Invoke(multiplier, new object[] { l, r });
            }
        }
    } 

その後、Nyaruruさんの詳しい記事を見つけたのでそれも投稿しておいた。Nyaruruさん(のさらに引用先)曰く、Expression Treeはセマンティクスをデータ構造として保持するからコードを構造を持った抽象的対象として扱えるのが利点で、作ってすぐさまコンパイルして動的メソッド生成のために使うなんてのは、邪道であるらしいのですが、まあCodeDOM使うとこんなに酷い事になりますよ、というお話でした。

ていうかこのコードはMult()呼ぶたびにコンパイルするのも酷いな。せめてstaticなコンストラクタでやるべきだ。

回ってフワフワするようにした

写真じゃ分かり難いけれど(当たり前だ)、z軸中心の回転とz軸方向の平行移動のアニメーションを加えた。さながらデジタルフィギュア。

コードはいつもの所においてある。実行ファイル入りなのでとりあえず遊びたい人もどうぞ。

---

追記：

ニコニコ動画にアップした。

体調

今までは疲労感＝倦怠感であって、明確で強い筋肉痛というのは週に二日くらいだったのだけど、いつの間にか完全に日常化してしまっている。腕から足から酷い。それに目がショボショボして目蓋が重く集中力が極端に阻害される。この数ヶ月で良くなったのは睡眠のリズムが整った事だけど、それでも毎日徹夜明け状態は変わらない。むしろ辛ければ寝ていれば良かった時期よりも一日の過ごし方に工夫が要るのでかえって耐え難い。だからミクさんでARというのも何ともアレだが、まだできる事があるというのは良い事だと思おう。

で、ミクさんをWPFのStoryboardでアニメーションさせるのが思ったよりも難航している。まあ気長にやるさ。