bashで3Dプログラミング

Posted on

今日はbashと基本的なシェルコマンドだけで3Dレンダリング(ワイヤフレーム)をしてみました!

Github: ledyba / Bash3D

bashプログラミングテクニック

今回つかったbashのプログラミングテクニックについて、記録のついでに少しまとめます。シェルスクリプトの基本的な所は書きません。

あと、3Dプログラミングに関しては本に書いてある事以上のことが書けそうにないので、(今は)書きません。

bashのデータ構造

bashでそこそこ本格的にプログラミングする事を考えた上でまっさきに「辛そう…」ってなるのはデータ構造ですよね!今回のソフトでは、ベクトルや行列を綺麗に扱えないと困ってしまうでしょう。まずはそこを見ます。

bashの「値」は基本的に全てただの文字列ですが、「リスト構造」にだけはネイティブで対応しています。こんな感じ。

$ lst=(1 2 3)
$ echo ${lst[0]}
1
$ echo ${lst[1]}
2
$ echo ${lst[2]}
3
$ echo ${lst[3]}
#空行

データにアクセスする時に、必ず${}でくくらないといけないのが少し面倒です。ちなみに添字を付けないとリスト全体…と思いたいところですが、最初の要素が出てきます。

$ echo $lst
1

全体を表示するには、${lst[@]}とします。

$ echo ${lst[@]}
1 2 3

リストを定義するときに使ったカッコが出てこないので、リストをコピーしたい場合はそのカッコを補わないといけません。

$ lst2=(${lst[@]})
$ echo ${lst2[@]}
1 2 3

慣れればなんとかなるレベルですが、面倒です。

さらに気を付けないといけないのは、”(1 2 3)”という文字列を返す関数を$(command)で受けて変数に代入してもリストになりません。

$ function fun() { echo \(1 2 3\); }
$ lst=$(fun)
$ echo $lst
(1 2 3)

ただの文字列として代入されます。どうしてもしたい場合は、evalを使わないといけません。

#文字列展開が起こってから実行されるので、lst=(1 2 3)というコマンドが実行されることになる
$ eval "lst=$(fun)"
$ echo ${lst[1]}
2
$ echo ${lst[@]}
1 2 3

通常の用途ならそれで済むのですが、今回は出来る限り高速に処理したいので、こうやってコマンドが沢山必要になるのは避けたいところです。

ところで。

後に言うダイナミックスコープとevalを組み合わせると連想配列みたいなことも出来るのですが、今回のプログラムでは基本的にこのリストでデータ構造を管理しています。

インデックスを指定したアクセスをするためにコストが高いevalを行わなくてよいこと、ベクトルと行列が主役なのでリストで十分な事が理由です。

$ matrix=(MAT 1 0 0 0 0 1 .....)
$ vector=(VEC 0 0 0 1)

最初の要素に型の名前を入れて適宜assertすることで、動的型検査みたいなことをしています。

変数のスコープ

名前空間は基本的にひとつしかありません。関数のネストは一切関係なく上書きされていきます。

$ val=0
$ function fun() { val=1; }
$ fun
$ echo $val
1 #関数内の代入だが、グローバル領域が上書きされてしまった
$ function fun2() { val2=2;echo $val2; }
$ fun2
2
$ echo $val2
2 #オフサイドルールではなく、存在しない変数はグローバルに作られるみたい

一種のダイナミックスコープみたいなものだと理解しておけば大体書けます。

呼び出し元の変数がうっかり書き換えられてしまう可能性があるのですごく不便ですが、逆に利用することも出来ます。

「関数の引数に変数名を渡すと、evalと組み合わせることで指定した変数に結果を格納する関数」が作れます。こんな感じ:

$function setOne() { eval "$1=1"; }
$ setOne val1
$ echo $val1
1
$ setOne val2
$ echo $val2
1

これを使って「最初の引数は入力、最後の引数は出力」と規約を決めて関数を設計しています。たとえば行列の計算はこちらにあります。見た目がアセンブラみたいになりますネ。

「こんなスコープ再帰が出来ないのでは?」という気がしますが、案外なんとかなります。というのも、$(command)はサブシェルとかいうので実行されるらしくて、さらにサブシェルは別の変数空間で動作するようです。つまり、このスコープの問題は起こりません。…たぶん。よく分からない(◞‸◟ )

 #再帰関数といえば階乗である
function kaijo() {
    arg=$1
    next=$(expr $arg - 1)
    if [ $arg -gt 1 ]; then
        result=`kaijo $next`
        echo $(expr $arg \* $result)
    else
        echo 1
    fi
}

kaijo 10 --> 3628800

ついでにevalとこのスコープの仕様を使うと、擬似的な連想配列を作れますが、値のコピーが難しいこともあって使いませんでした。

$ val_key1=123
$ val_key2=456
$ function showKey1() { echo $(eval "echo \$$1_key1"); };
$ showKey1 val
123

固定小数点とexprコマンド

シェルスクリプトで計算をするには、exprコマンドを使う必要があります。このコマンドは整数の計算にしか対応していません。

$ expr 1 \* 2 \* 3
6
$ expr 1 \* 2 \* 3.2
expr: non-numeric argument

3Dソフトでそれでは困るので、整数だけで小数を計算します。それにはどうするかというと、扱う小数を適当に1000倍とかして、それを小数に見立てます。つまり、1000は1、1001は1.001、3140は3.140だと思うことにします。その場合、計算する時にちょっとした手間が必要になります。

#3.14 + 3.14
$ expr 3140 + 3140 #足し算・引き算はそのまま足し算で大丈夫です。
6280
#3.14 / 3.14
$ expr 3140 \* 1000 / 3140 #割り算するときは、倍率で掛けてから割ります。この世界では「1.0」は1000であることに注意
1000
#3.14 * 3.14
$ expr 3140 \* 3140 / 1000 #掛け算するときは、掛け算してから倍率で割ります。
9859

倍率を1000じゃなくてもっと高くすれば、もっと小数演算の結果が正確になります。

私の環境だけかもしれませんが、exprはどんな大きな数でもちゃんと扱ってくれるようなので、ガンガン倍率を上げていけば誤差がまったく気にならなくなります。これはすごい楽…(楽すぎた)

なお、sin/cosは(テイラー展開を自分で展開しない限り)計算できないようだったので、0~360度について事前に求めておいて、これを先ほどのリスト構造にいれて使っています。平方根も使いたかったんですが、どうやっても無理そうなので諦めました。

表示:オフスクリーンレンダリング

もちろんターミナルにはピクセルがないので、ターミナル上に表示する文字で代用します。普通に表示するなら、tput cup (line) (col)を使うとターミナル上の任意の位置に移動できるので、これを使いたいところです。

# 斜め45度に線を引いてみた
for1;do
    for2;do
        tput cup y x
        echo -n \*
    done
done

が、この方法では線を引く様子が目で分かるぐらい時間がかかります。これで最初レンダリングしたみたのですが、レンダリング中に線を引いている様子がバッチリ見えてしまい、あんまり回転アニメーションしてるように見えませんでした。

ほんもののゲームでもそういう、更新途中の絵が表示されてしまってちらつく問題があります。ほんもののゲームでは、そういう時は、書いてもすぐには表示されない「裏の画面」を用意しておいて、そちらに時間を掛けてレンダリングして、最後に表示される「表の画面」と入れ替える「オフスクリーンレンダリング」と「ダブルバッファリング」をしているのでした。私も同様にそのテクニックを使いました。

具体的には、こんな感じ。

SCR_BUFF=$(head -c $size < /dev/zero | tr '\0' '_')

今回の「ピクセル」は文字ですから、「裏の画面」はただの文字列になります。
画面のサイズを取得しておいて、そのサイズ分何もない文字列を意味する「_」をtrを使って入れておきます。空白でいいじゃないかと思うかもしれませんが、シェルスクリプトでは空白はかなり重要な意味を持っていて面倒なのでこうしています。

「ドットを打つ」時は、対応する文字列中のインデックスの文字を*に書き換えます。

SCR_BUFF=$(echo $SCR_BUFF | sed s/./*/<置き換える文字の位置>)

zshだとstr[12]=\*ってやると文字列を更新できるのですが、bashでは使えなかったので、sedを使って書き換えています。この処理が非常に重く、今回のレンダリングエンジン(?)の処理時間の殆どはこのドットを打つ処理に費やされています。上手い方法を知ってたら教えてください!

最後に表の画面に更新する方法ですが、単純にechoするだけです。この時、何も入っていない所には_が入っていたので、空白に直します。

$ tput cup 0 0
$ echo -n $SCR_BUFF | tr '_' ' '

ちゃんと幅に合わせて改行されるので、うまく表示できます。

高速化

インライン展開 → プロセスforkを抑える

当初綺麗に関数を使って書いた所、行列の掛け算に1秒近く掛かってしまい、さすがにレンダリングは出来ないのかな…と諦めかけたのですが、行列やベクトルの掛け算や足し算を全てインライン展開したところ、座標計算に関しては十二分すぎる程の速度が出せるようになりました。前述した通り、今回遅いのはオフスクリーンにピクセル(文字)を書き込む処理が異常に時間が掛かっているからです。

やはりプロセスの起動はものすごい重い処理のようで、可読性をそこまで損なわない範囲でプロセスの呼び出しを一箇所にまとめています。

感想

線を引くのにえげつない時間が掛かるのを見た時、初めて触ったBASICの処理系を思い出しました。懐かしい。

ところで、なんで3Dライブラリを作るとOpenGL風の名前にしちゃうんだろ。

ニコ動に「次はシェーディングだな!」というコメントがついているのですが、テクスチャの表示をやってみたいです!7色くらいしか使えないので、たいへんそう…。影は「暗い色」が基本的に無いのでかなりつらい気がします。

  1. x=0;x<10;++x []
  2. y=0;y<10;++y []

3Dをゼロからレンダリングするテスト

Posted on

C++から動画を直接出力する自作ライブラリ「Dolly」の上で、フルスクラッチで3Dのレンダリングを行なっています。

とりあえずは、ワイヤーフレームの正方形をぐるぐる。

テクスチャもなんとか実装しましたが、三角形の継ぎ目が出ているのが気になります。でもどこが悪いのかわからん…

このサンプルのソースコードはこちらです:

https://github.com/ledyba/Dolly/blob/master/src/sample/super_px.cpp

東北地方太平洋沖地震の地震波を音にしてみた

Posted on

地球物理をやる学科に居るので、地震も学んでいるのですが、ふと思いついたので地震計の観測した波形を音に見立ててWAVファイルにしてみました。

気象庁の強震波形(平成23年(2011年)東北地方太平洋沖地震)のページから、福島県郡山市朝日の、震度6のデータです。なんだか、雷みたいな音がしますね!

ソースコードはこちらです。PythonでCSV読んでWAVに書いているだけですが、一応使い方のサンプルとしては使えるのではないでしょうか。

凪のあすからの恋愛関係からページランクを求めてみました

Posted on

最近友人に薦められて、今日が最終回らしい「凪のあすから」を見ています。今6話なのですが、恋愛関係が大変複雑です。「深夜の昼ドラ? アニメ『凪のあすから』の人物相関図つくったよ」というページがあったので、引用しますと、こんな感じ。

うーん、グラフ理論って感じだ。そう思ってからピクシブの非公式カップリングのページを見てしまうと、もはや隣接行列か何かにしか見えなくなってきた…。

20130403隣接行列…グラフ理論…はっ!ページランク求めようページランク!

というわけで、今まで一回ぐらいやってみようと思っていた、ページランクの計算を実際にやってみました。こちらの「Google の秘密 – PageRank 徹底解説」を参考にしました。人間をページに見立てて、好意をリンクに見立てます。

ページランクとは結局何なのか

昔私が読んだ本だと、

  • 全てのページはページランク値PRを持っていて、
  • 他のリンクしているページそれぞれににPR/N(Nはページ内のリンク総数)を与えており、
  • ページランクの値PRはリンクされた他のページから貰ったページランク値の和に等しい

というのを読んだ気がします。つまり、貰ったページランク=あげたページランク=そのページのページランクという等式が成り立つということです。

元の論文にもそんな図が載っています。

20130403-2

Lawrence Page, Sergey Brin, Rajeev Motwani, Terry Winograd, ‘The PageRank Citation Ranking: Bringing Order to the Web’, 1998,
http://www-db.stanford.edu/~backrub/pageranksub.ps

これをグラフ理論と確率過程の話で表現すると、こうなります。

  • 各ページから他のページへ移動する確率の表を用意する(とりあえずすべてのリンクで等確率=1/リンク総数とする)
  • その表に従ってあるページから出発してリンクをたどり続ける
  • すると、どのページから出発しても、N回リンクをたどった時にあるページPに居る確率がN+1回目のそれと変わらなくなってくる(これを定常状態と呼びます)
  • その時の、N回目でそれぞれのページに居る確率がページランクです。

これがホントに元の話と一緒なのか?ということなのですが、

  • “もらうページランク”=(自分も含めた他のページにN−1回目に居る確率*それぞれの遷移確率)の和=N回目にそのページにいる確率
  • “あげるページランク”=(N回目にそのページにいる確率 * そのページから他のページへの遷移確率)の和=N回目にそのページにいる確率1

なので、もらうページランク=あげるページランクという関係は確かに成り立ちます。

求めてみる

ページ遷移確率表=行列

実際にもとめてみましょう。まずは、鍵となるページ遷移確率表を用意します。さっきの表を見ながら、ひーくんからまなかへ遷移する確率は1…みたいな感じで行列を作るとこうなります。一応、今回は友情は無視して恋愛関係だけ考えます。

M =
 ひーくん ちさき  まなか   要    さゆ   美海   つむぐ  あかり  至    みをり ←好き ↓好かれる
   0.00000   1.00000   0.00000   0.00000   0.00000   1.00000   0.10000   0.00000   0.00000   0.00000 ひーくん
   0.00000   0.00000   0.00000   1.00000   0.00000   0.00000   0.10000   0.00000   0.00000   0.00000 ちさき
   1.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.10000   0.00000   0.00000   0.00000 まなか
   0.00000   0.00000   0.00000   0.00000   1.00000   0.00000   0.10000   0.00000   0.00000   0.00000 要
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.10000   0.00000   0.00000   0.00000 さゆ
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.10000   0.00000   0.00000   0.00000 美海
   0.00000   0.00000   1.00000   0.00000   0.00000   0.00000   0.10000   0.00000   0.00000   0.00000 つむぐ
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.10000   0.00000   0.50000   0.00000 あかり
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.10000   1.00000   0.00000   1.00000 至
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.10000   0.00000   0.50000   0.00000 みをり

普通、ページには複数のリンクがありますが、このアニメだと思い人が二人以上居るのは不倫と勘違いされて殴られた至だけです。

あと6話時点ではつむぐは誰も好きではないので、つむぐの行が全部0になるはずなのですが、そうしてしまうとページランクが求まりません2

これを元論文ではdangling pageと呼んでいます。元論文では「あんまり重要じゃないから削除して計算しています」と無視してしまっているのですが、つむぐくんを無視するのはかわいそう(?)なのと、仮につむぐくんを消してしまうと今度はまなかがdangling pageになって…となって最終的に大人組しか居なくなってしまうので、次に紹介するrandom surferモデルに習って、みんなを平等に好きなことにしました。

計算は掛け算するだけ

元の定義通り、「あるページからスタートしてリンクをたどり続けた時に、それぞれのページに居る確率が一定になってくる」まで遷移させ続けます。

たとえば、ひーくんからはじめましょう。はじめはひーくんに居ると仮定しているので、ひーくんに居る確率は1で、他みんな確率0とします。

1回リンクをたどると、グラフを見るとひーくんはまなかが好きなので、まなかの確率が1になるはず。

これを我々が先ほど作った行列を使って計算しようとするなら、先ほどの行列Mに、ひーくんを1としたベクトルPを掛ければよいです。

octaveでやるとこんな感じです。

octave:95> M=[ 0 1 0 0 0 1 1/10 0 0 0 ; 0 0 0 1 0 0 1/10 0 0 0; 1 0 0 0 0 0 1/10 0 0 0; 0 0 0 0 1 0 1/10 0 0 0 ; 0 0 0 0 0 0 1/10 0 0 0 ; 0 0 0 0 0 0 1/10 0 0 0; 0 0 1 0 0 0 1/10 0 0 0 ; 0 0 0 0 0 0 1/10 0 1/2 0; 0 0 0 0 0 0 1/10 1 0 1; 0 0 0 0 0 0 1/10 0 1/2 0 ]
(略)
octave:96> P = [1; 0; 0; 0; 0; 0; 0; 0; 0; 0] ←ひーくんに居る確率を1としたベクトル
octave:97> M*P
ans =

   0
   0
   1
   0
   0
   0
   0
   0
   0
   0

というわけで、まなかは三番目の要素で、それが1となったので、見事正しい結果を得られています。これを沢山続けてみましょう。10000回ぐらい掛けると変わらなくなってきます。

octave:98> (M^10000)*P
ans =

   0.00000
   0.00000
   0.00000
   0.00000
   0.00000
   0.00000
   0.00000
   0.22222
   0.55556
   0.22222

…あれぇ?最後の3人、つまり大人組以外ゼロ?なんか間違ったかなぁ?

ランダム・サーファー・モデル

こうなってしまった原因は、大人組から子供組へ繋がっていないことです。つむぐくんは全員が好きだと仮定したので、子供組から大人組へはつむぐくんから遷移することが出来るのですが、その後は大人組の間でぐるぐるぐるぐる巡り続けることになります。ひーくんからはじめても長い目で見ればいつかは大人組に入ってしまい、そこで定常状態となるので、子供組の確率がすべて0になってしまったのです。グラフ理論の言葉でいうと、グラフを辿ると元の位置に必ず戻ってこれる大人組のような部分を再帰類と呼び、そうではなく元に戻れない子供組のような部分を非再帰と呼ぶそうですが、あまりこれ自信をもって解説できないので、しません…。

元の論文ではこれをRank sinkと読んでいて、そのためにランダム・サーファー・モデルというのを考慮しています。これはどういうモデルかというと、「時々はまったくリンクされていない別のページに飛ぶ」と考慮しているモデルです。具体的には、新しい遷移確率行列を次のように定義します。

octave:99> M2=(M*1/1.15)+((1/10)*0.15/1.15)

これは、まず今までの確率を少し下げた上で、あるページから他のページへ僅かな確率で飛ぶようにした、ということです。凪のあすからの言葉を使えば、ひょっとするとひーくんは他の女の子(や男の子)にも興味を持つかもしれないよねというモデルです。

この修正したモデルでは、あるページから他のページへ少ない確率ながら必ずジャンプできるので、Rank sinkのような問題は起こりません。

早速、計算してみましょう。

octave:100> (M2^10000)*P
   0.095976 ひーくん
   0.060684 ちさき
   0.106569 まなか
   0.043208 かなめ
   0.023111 さゆ
   0.023111 美海
   0.115780 つむぐ
   0.135981 あかり
   0.259599 至
   0.135981 みをり

やはり両思いは強いのですね!(?) 子供組を見ると、やはり特に重要そうなキャラクターであるひーくん、まなか、つむぐの三人のページランクが高いので、おおよそページランクの妥当性を示せたと言っていいんでしょうか!?とりあえず、この三人に注目して今後は視聴を続けたいと思います!

蛇足:固有値・固有ベクトルを使って求める

なんか10000回も行列を掛け算するのは…という感じしますか?元の論文だとその方法で計算してると書いてあるのですが、これは行列の固有値・固有ベクトルを使っても求めることが出来ます。

固有値・固有ベクトルをおさらいしておくと、

  • ある行列Mに対する、
  • 固有値μ・固有ベクトルVは、
  • 次のような関係を満たすものです
  • MV = μV

ここで、固有値μ=1なら、MV=Vとなって、N回目とN+1回目でのそれぞれのページに居る確率が全く同じという事になり、ちょうど私達が求めたかったベクトルと一致しますね!

octaveで求めるとこんな感じになります。

octave:142> eig(M2)
ans =
   1.00000 + 0.00000i ←求めたい、固有値1の固有ベクトルは1番目だと分かる
   0.78950 + 0.00000i
   0.19713 + 0.55213i
   0.19713 - 0.55213i
  -0.86957 + 0.00000i
  -0.59237 + 0.00000i
  -0.25222 + 0.45314i
  -0.25222 - 0.45314i
   0.00000 + 0.00000i
   0.00000 + 0.00000i

octave:143> [V, D] = eig(M2); ←固有値と固有ベクトルを再度格納
octave:144> V(:,1) ←固有ベクトル1番目を取得する
ans =
   0.252078
   0.159383
   0.279899
   0.113484
   0.060701
   0.060701
   0.304091
   0.357146
   0.681825
   0.357146
↑固有ベクトルは整数倍しても固有ベクトルで、octaveはベクトルの大きさを1に正規化するので、さっきの結果と合わない
octave:146> V(:,1) ./ sum(V(:,1)) ←ベクトルの要素の合計が1となるようにすると…
ans =
   0.095976
   0.060684
   0.106569
   0.043208
   0.023111
   0.023111
   0.115780
   0.135981
   0.259599
   0.135981
↑先ほどの結果と一致する!

こんな感じで計算できます。固有値ってこんなのにも使えるんですね〜。ちなみに、この行列の固有値のうち、絶対値が一番大きいものはかならず1になることがペロン=フロベニウスの定理から言えるらしいのですが、よくわからないので解説できません…。

感想

凪のあすからは線形代数と確率過程とグラフ理論を学べる神アニメ

ていうか実は前期のアニメで一番好きだったのは「いなり、こんこん、恋いろは」で、京都と出雲にまで行ってきたのですが、それを書く日は来るのだろうか…!?

  1. 他のページへの遷移確率の和は1です。リンク総数*(1/リンク総数)=1。 []
  2. つむぐから他のページに遷移できないので、確率がそこで皆0になってしまう []