球オブジェクトが回転しても天地が回転しないテクスチャをつくる

はじめに

球体POVというものを作成しています．

tajmahal0707.hatenablog.com

こんな感じで球体に映像を表示することができるガジェットで

このように完全にボールの中に収納してバッテリー動作することができます．

球体POV

このようなガジェットを制作してると，当然やりたくなることがあります．

そう

転がしても天地が同じ（画像は回転しない）になるボール

近いイメージはこんな感じ？

映画『ジュラシック・ワールド』ボーナスDVD特典映像

縁日のおもちゃとかでボールの中に液体があり，中の物が浮いているので転がってもその向きが変わらないおもちゃとかありませんでした？　名前が理解らなくて検索できませんでした．

やってみたいと思いませんか？　　私はそう思います．．

そこでそんなこともあろうかと，このガジェットのマイコン部分には

M5Stack ESP32 Mpu9250 9Axies Motion Sensor Core Development Kit Extensible IoT Development Board for Arduino

出版社/メーカー: M5Stack
メディア: エレクトロニクス
この商品を含むブログを見る

を採用しています．

内部にMPU9250を搭載し，センサフュージョンから自己の姿勢quaternionやオイラー角を算出することができます．

これを使ってボールが回転しても，表示する映像は回転しない機能を実装しようと思います．

球体テクスチャの基本

この球体POVは，画像でテクスチャを与えるとそのLEDの位置に相当する色をテクスチャから取得して画像を表示します．

そのときのテクスチャを扱いやすい形式で与えることで正しい像を表示することができます．ここでは貼り付ける画像を「正距円筒図法」に基づいたものと定義しています．

ja.wikipedia.org

例えばこんな画像です．

f:id:tajmahal0707:20191126111738p:plain
これは，パノラマ画像とも呼ばれ，近年のCG分野で球体に貼るテクスチャとしては比較的一般的なもので、メルカトル図法とは縦の比が異なる親戚のような図法です．（地図の表現法について詳しくはこのへんあたりを参照）

この画像は基本的にはよくネットでも拾ってくることもできますし，thetaなどの３６０度カメラでも出力モードに設定されているので，入手することはそれほど苦労することは無いかと思います．

この図法の特徴は

緯度と経度がそのまま画像上の座標値（uv座標）になっている
緯度と経度で作る長方形が正方形になる

一方で北極・南極に向かうほど歪みが大きくなってしまうという欠点があります．

この画像テクスチャを球体にマッピングすることを考えます．

通常はこの処理はGPUで行ってしまうのですが，マイコンで処理を行うために一つ一つばらばらにしてみます．

上の地球の図面を例に考えると，テクスチャ座標を（u, v）（ここで画像サイズに依存しないように正規化して｛0.0 ≦ u, v ≦ 1.0｝とする）とします．

つまり，左上が（u, v）=（0.0, 0.0）で，右下が（u, v）=（1.0, 1.0）となります．（ただし実画像では緯線（±90度）と経線方向（360度）の長さ比は1:2で経線方向が長い）

赤道が北緯0度．北極は北緯９０度（ $-{\pi}/2$ ），南極が南緯９０度（ ${\pi}/2$ ）となっています．東経は360度で地球一周（ $2{\pi}$ ）することとなり，u軸方向の画像中央は東経0度（グリニッジ）上で、画像中央（東経0度、北緯0度。アフリカの海のど真ん中）のuv座標は（0.5, 0.5）となる．

例えば東京を例（北緯36度，東経140度）に取ると

$u=140/360+0.5=0.888$ 角度成分で正規化して，グリニッジ分（0.5）を足す
$v=0.5-36/180=0.3$ 角度成分で正規化して，赤道分（0.5）から引く

となり，東京のuv座標(u, v)＝(0.888, 0.3)は

f:id:tajmahal0707:20191118145553p:plain

だいたいこの辺になるということです．

このように球面上の座標が判れば，換算したuv座標から対応する座標の色を算出できることになります。

球体POVの画像生成

さて、球面のテクスチャ座標が分かったところで、球体POVのLED各色を算出する方法について考えます。

球体POVは

球体POVデモ

最初の10秒くらいをみると分かるように円弧状に配置したLEDを回転させ、角度ごとに表示する色パターンを変えることで残像効果（Persistence of vision）によって映像を作る方法です。

そのため回転軸を北極南極と見なすと先のパノラマ画像のuv座標と相性の良い方法（LED配列方向がv軸方向、回転方向がu軸方向に対応）なのです。なぜならテクスチャから縦方向に順に色を読んで配置していくだけでLEDの色を決定することができるのですから。

　　もちろん転がらない場合には．．．ね

この球体POVはLED配列方向に45個のLEDを、回転軸方向を120分割することによって120×45の分解能を持つディスプレイになります。

地球を表示するとこんな感じ（まだLEDの問題でフルカラーで出せていませんが）

球体POV

これを物体の回転に合わせてテクスチャ座標も回転させてみようというのが今回の試みです。ここからはまずprocessingを使ったシミュレーションで。

テクスチャ回転①

シンプルにテクスチャを横回転させてみましょう。

北極南極の軸は固定のままで、極軸周りに回転。オブジェクトは不動です．

longitude
スタートとなるuの座標が変化するだけなので簡単ですね。

f:id:tajmahal0707:20191119075538p:plain

青矢印のようにLEDスキャンを開始する位置をずらすことで，オブジェクトは動いていないのに映像が回転しているように見えることがわかります．

テクスチャ回転②

今度は縦回転。ｘ軸周りに回転するので北極南極軸が回転します。

何も考えずに開始するvの座標値を変えてみましょう。

latitude

あれあれ、まあそうなりますよねえ。。。

極点周りで座標系が歪んでしまっているのにそれを考慮せずに回転させたらこうなります。

そこで天頂位置を変動させる方法を考えます。最初、この計算について１週間以上悩みました。。最初メルカトル図法のグーデルマン関数あたりの方に向かってしまい、相当ハマりました。

テクスチャ回転③

ポイントは「大円」です。

ja.m.wikipedia.org

大円とは球をその原点を通る平面で切った時にできる3次元円のことでその半径は球の半径に等しくなります．

f:id:tajmahal0707:20191119111806p:plain
球体POVの映像生成方法は，天頂軸の北極・南極を通る大円（以後LED大円，図中緑の円）に北極から南極に向けた半円上にLEDのを配置します．これがテクスチャuv座標系におけるv軸方向の座標を決定します．

u軸方向の座標は，LED大円を天頂軸回りの大円（以後回転軸大円，図中の青い円）に沿って回転させることでu軸をスキャンします．図からも分かる通り，回転軸大円とLED大円の法線ベクトルは直交します．

そのときの各LEDの座標値をuv座標に変換することで回転するテクスチャ値を得ることができると考えられます。

ここでもう一つのポイントはquaternion（四元数）です。

quaternionとはある回転軸ベクトル ${v}$ 周りに $\theta$ だけ回転させた量で回転を記述する方法です．

まずベクトルvを ${v'}=\dfrac{v}{\sqrt{v_x^ 2+v_y^ 2+v_z^ 2}}$ と単位ベクトル化します．

回転角をθ(rad)としたとき，この回転を表すquaternion：qは次式で表されます．

${q}=(q0, qx, qy, qz)=(cos(\dfrac{\theta}{2}),　{v'} sin(\dfrac{\theta}{2}))$

これをLED大円および回転軸大円の回転軸に対して適用します．

LED大円におけるLEDの配置は，北極から南極まで $-{\pi}/2〜{\pi}/2$ の180度を45分割します．この座標値はquaternionを用いると簡単に記述できます．

これを大円周りの回転に適用すると，原点を通るLED大円の法線ベクトルを ${n}_{LED}$ とすると，この法線ベクトル ${n}_{LED}$ がquaternionの回転軸ベクトル成分となります．そして北極から南極に向かって $180/45=4$ 度間隔で回転させていくと４５個のLEDの各座標値を得ることができるのです。これによってある経度におけるLED座標を求めることができ、それを天頂軸で回転させます。

天頂軸回転もquaternionで簡単に記述でき、天頂ベクトル ${n}_{apex}$ を回転軸ベクトルにした回転を $180/120=1.5$ 度間隔で回転させていくことに相当します。

天頂ベクトルのセンサによる回転

上記までである天頂ベクトルを与えた場合のLEDのuv座標を計算する準備ができました．

回転がない場合の天頂ベクトル ${n}_{apex}$ は単位球のy軸上の天頂（北極に相当）に位置し，その座標値は ${n}_{apex}=(0, 1, 0)$ となります．またLED大円の法線ベクトル ${n}_{LED}$ は ${n}_{apex}$ と直行するx軸上と考えることができ，その座標値は ${n}_{LED}=(1, 0, 0)$ とします．

そこでMPU9250で取得した姿勢情報から天頂ベクトル ${n}_{apex}$ を回転させることを考えます．

MPU9250でquaternionを求める方法はいくつかありますが，M5stackのサンプル例にある「MPU9250BasicAHRS」を用いるのが一番簡単です．

f:id:tajmahal0707:20191119123212p:plain

このサンプルの322~325行目にquaternionを取得する部分があります．

f:id:tajmahal0707:20191119123626p:plain

オイラー角（roll, pitch, yaw）も求められますが今回は使いません．

センサから取得したquatrenion を $q_{sensor}=(q0,qx,qy,qz)$ とします．

$(qx, qy, qz)$ が回転軸ベクトル． ${q0}$ が回転量に相当し，回転軸ベクトル $(qx, qy, qz)$ は単位ベクトルであるとします．

このquaternion $q_{sensor}$ にもとづいて 2つのベクトル（ ${n}_{apex}$ ， ${n}_{LED}$ ）を回転させます．この2つのベクトルはそれぞれ軸方向の単位ベクトルなので0の項が多く計算を簡単にできそうですね．

通常のベクトル ${v}=(v_x, v_y, v_z)$ をquaternion ${q}_{sensor}=(q0, qx, qy, qz)$ で回転させる回転行列 ${R}_{sensor}$ を求めると

${R}_{sensor} = \begin{pmatrix} 1-2(q_y^ 2 + q_z^ 2)&2(q_x q_y - q_z q_0)&2(q_x q_z + q_y q_0) \\ 2(q_x q_y + q_z q_0)&1-2(q_x^ 2 + q_z^ 2)&2(q_y q_z - q_x q_0) \\ 2(q_x q_z - q_y q_0)&2(q_y q_z + q_x q_0)&1-2(q_x^ 2 + q_y^ 2) \end{pmatrix}$

天頂ベクトル ${n}_{apex}$ を回転 ${R}_{sensor}$ で移動させると

$\begin{eqnarray}{n'}_{apex} &=& {R}_{sensor}\begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix} \\ &=& (2(q_x q_y - q_z q_0),　1-2(q_x^ 2 + q_z^ 2),　2(q_y q_z + q_x q_0))\end{eqnarray}$

と簡単な式になります．零があるのでたくさん消えますね．

例えば天頂ベクトルをLED大円の法線ベクトル ${n}_{LED}$ で90度（ ${\pi}/2$ ）回転させることを考えると，このときのquaternion ${q}_{LED} = (cos(\pi/4), sin(\pi/4), 0, 0)$ から天頂ベクトル ${n}_{apex}$ を回転させると，

${R}_{LED} = \begin{pmatrix} 1&0&0 \\ 0&1-2(sin(\pi/4)^ 2)&2( -cos(\pi/4) sin(\pi/4)) \\ 0&2 cos(\pi/4) sin(\pi/4)&1-2(sin(\pi/4)^ 2) \end{pmatrix}$

$\begin{eqnarray}{n'}_{apex} &=& {R}_{LED}\begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix} \\ &=& (0,　1-2(sin(\pi/4)^ 2),　2 cos(\pi/4) sin(\pi/4)) \\ &=& (0,　0,　1)\end{eqnarray}$

となり，天頂ベクトル ${n}_{apex} = (0,1,0)$ が90度回転して $(0,0,1)$ へと変換され，x軸回りにy軸を90度回転させてz軸になったことがわかります．

同様にLED大円の法線ベクトル ${n}_{LED}=(1,0,0)$ を回転 ${R}_{sensor}$ で移動させると

$\begin{eqnarray}{n'}_{LED} &=& {R}_{sensor}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} \\ &=& (1-2(q_y^ 2 + q_z^ 2),　2(q_x q_y + q_z q_0),　2(q_x q_z - q_y q_0))\end{eqnarray}$

これも簡単な式にすることができました．

それぞれ計算が簡単なので，別々の関数にしたほうが良さそうですね．

また，Quaternionには球面線形補間という性質があり，2つのベクトルを円弧状に補間することができます．

LED大円上の点p0とp1があるとき，パラメータt（0.0 ≦ t ≦ 1.0）を用いてその間を補間する中間のベクトルp(t)を求めることができ，その式は

$p(t) = \dfrac{sin((1-t)\theta)}{sin(\theta)} p_0 + \dfrac{sin(t\theta)}{sin(\theta)} p_1$

と表すことができます．

しかし，ここではこのあともセンサによる回転などが加わるのと，実装時にはこのあたりの計算はすべてすっ飛ばしてしまうことになるので，いろいろ検討したのですが今回は採用しないことにしました．

球上の座標の緯度経度計算

上記までセンサによって回転した天頂ベクトル ${n'}_{apex}$ を求めることができました．この天頂ベクトルは単位ベクトルなのでそのまま球上の座標値として使うことができ，この座標値をもとに回転した天頂の緯度経度 ${n'}_{apex}$ （ひいてはテクスチャ上のuv座標）を算出します．

直交座標から極座標への変換公式を用います．単位球なので $r=1$ となり，

$\begin{eqnarray}\left\{\begin{array}{l}r &=& 1 \\ \theta &=& arccos({z}) \\ \phi &=& arccos(\dfrac{x}{\sqrt{x^2+y^2}})\end{array}\right.\end{eqnarray}$

で求めることができます．

この極座標変換の結果から緯度： $\theta$ と経度： $\phi$ が算出できます．この $\theta$ と $\phi$ は単位がRadianなので，PIで割って正規化することでuv座標に変換することができるのです．

センサからオブジェクトの回転を踏まえたLEDの各点の球面上の座標を求め，その座標値の緯度経度ーひいてはテクスチャのuv座標を求めることで，各LEDのカラー値を求めることができます．

そこで，m5stackから取得したセンサ値をシリアル通信でprocessingに送り，その値をもとにテクスチャ値（というかLEDの位置にcubeを配置してそのcubeの色を指定）をシミュレーションしてみました．

Rotation test

画面の下がセンサ値に応じて回転している球体（正確には球表面上にcubeを配置したもの）で，その各LEDに相当するcubeの色を上のパノラマ画像のテクスチャ座標から得ています．（センサ値と実際の回転が一致していないのはご愛嬌．）

そして上のパノラマ画像上には，テクスチャ値を取得するために使用したuv座標の等緯線をマッピングしています．

動画を見るとわかりますが，オブジェクトの北極・南極が回転していろいろな位置に移動していることがわかります．そこを中心に新たにuv座標を取り直すことでパノラマ画像から回転に応じたテクスチャを作り直すことができるのです．

この処理の擬似コードはこちら

  // センサ値からquaternionを取得
  quat[4] = getQuaternionFromSensor();
  
  // 天頂軸回りに回転：モーターの回転方向
  for(int u = 0; u <= LONGITUDE; u ++){
        float uu = float(u) / float(LONGITUDE);

        // LED大円周りに回転：各LEDの色を求める
        for(int v = 0; v <= LATITUDE; v ++){
            float vv = float(v) / float(LATITUDE);

            // LED大円回転（vv*PI [rad]），結果をvec2に格納
            float[] vec2 = rotate(axisX, vv*PI, axisY);

            // vec2を回転軸大円で回転（uu*PI [rad]），結果をvec1に格納
            float[] vec1 = rotate(vec2, uu*PI*2.0 , axisX);

            // センサからの回転 quat でvec1を回転し，結果をvec0に格納
            float[] vec0 = rotateQuaternion(quat, vec1);
          
            float[] uv = pos2uv(vec0);            // vec0 のuv座標を変数uvに格納
            // uv 座標からテクスチャのカラーを取得
            color c = getColor(int(uv[0] * image_width), int(uv[1] * image_height));
        }        
    }

quaternionは回転させる順番で結果が異なるので，順番に注意．

次回はこれを実際に球体に落としこむ作業に入ります．

大量の掛け算，三角関数計算があるので，ダイエットがんばらないと．．