2016-06-01

PRMLを読みすすめるための記号一般を整理する

PRMLを読んでいるとギリシャ文字がよく現れます。知識の基盤がある人なら「ああ、アレね」と暗黙のパラメータのように読みすすめられるのかもしれませんが、私のレベルでは「以前も調べた気がするけどなんだっけコレ」となってしまいます。数式中に現れるギリシャ文字は慣例的で、例えばθは角度でiは虚数を意味するなどがありますが、馴染みがない記号が現れるとそれだけで躊躇してしまうことが多々あります。特に書籍を読み始めた当初はなおのこと。

よく使うものは計算ノートの一番後ろのページにメモ書きで残しているのですが、そろそろノートが一杯になってしまうので捨てた後でも振り返れるようにブログに残してみようと思います。色々と自分で調べて納得したものもありますが、私の場合は老化などの諸々な事情もあって物覚えが悪いというか、忘れっぽくなってきているようなので「なんだっけこれ？調べた気がする」ってときに振り返れるしおりが必要です。

というわけで今回のブログ記事はPRMLの内容そのものではなく、読み進めるための前提条件として必要になった事前知識の整理といったところでしょうか。幾つかの書籍をあたってみて学問分野ごとにできる限り一般的っぽいものを選んでみたつもりです。読み方はネットで調べてみましたが、ちょっと自信ありません。分野や現場によって色々あるらしいので。

統計学に関するもの

式	読み	意味
${\sigma}$	シグマ	標準偏差。標準偏差の自乗が分散なので、分散は ${\sigma ^ 2}$ として表されます
${\rho}$	ロー	回帰分析における相関係数
${\mu}$	ミュー	期待値
${\eta}$	エータ	分布の位置母数パラメータなど
${\lambda}$	ラムダ	分布の尺度母数パラメータなど
${\phi}$	ファイ	線形基底関数モデルにおける基底関数など
${\mit \Sigma }$	シグマ	分散共分散行列。和の記号と違って、この場合は斜体で表されます
${\propto}$	プロポーション	比例

日本で比例記号をプロポーションって読み上げる人っているのかちょっと謎ですね。

「(事後分布) ${\propto}$ (尤度)(事前分布)」って書かれていても「事後分布は事前分布と尤度の積に比例する」って読んだほうが早そうですし、実際に書籍を読んでる最中もそう読んでいます。

ベクトル解析に関するもの

式	読み	意味
${\delta}$	ラウンドディー	偏微分の記号
${\nabla f}$	ナブラ	fの勾配
${\nabla \cdot f}$		fの発散
${\nabla \times f}$		fの回転
${\nabla ^ 2}$		ラプラシアン
${\|\|V\|\|}$		ベクトルVのノルム。直観的に解釈するならVの大きさ

書籍では ${\nabla f}$ はスカラー場の勾配を簡単に表記するために使われます。ラプラシアンは線形基底回帰モデルでラプラス基底を使った例などで取り上げられます。発散や回転は書籍内では現れませんが、ベクトル解析の分野では基礎となるものらしいので、ここでまとめて列挙しておきました。

一般数学

式	読み	意味
${\lambda}$	ラムダ	行列の固有値など
${\ln}$	エルエヌ	自然対数(ネイピア数eを底とする対数)
${M ^ T}$	転置	Mで表される行列またはベクトルの転置
${M ^ {-1}}$	逆行列	行列Mの逆行列
${\Sigma}$	シグマ	和の記号
${\theta}$	シータ	角度。普通はラジアン単位
${\Pi}$	パイ	積の記号
${\tau}$	タウ	合成積の区間

${\tau}$ はこの分類でいいのかちょっと謎ですが、とりあえず忘れないように入れておきました。合成積というと馴染みありませんが、畳込みとかコンボリューションと言い換えればソフト屋さんなら通じやすいかと思います。 ${\Sigma}$ は今さら感があるのですが分散共分散行列と混同して困ってしまったので、区別する目的で記載しました。

確率と分布

式	意味
${P(X)}$	確率変数Xの分布
${P(X=1)}$	確率変数Xが標本値1をとる確率
${P(X,Y)}$	X,Yの同時確率分布
${P(X\|Y)}$	Yが与えられたもとでのXの条件付き確率
${E[X] }$	確率変数Xの期待値
${V[X] }$	確率変数Xの分散
${Bi(n,p)}$	二項分布。n=1の場合はベルヌーイ分布
${N(\mu,\sigma ^ 2)}$	期待値 ${\mu}$ と分散 ${\sigma ^ 2}$ を分布のパラメータとする正規分布。別名ガウス分布。
${Cov(X,Y)}$	確率変数XとYの共分散
$\displaystyle {P(X)=\sum^{}_{Y} P(X,Y) }$	加法定理
${P(X,Y)=P(Y\|X)P(X)}$	乗法定理
${P(B\|A)=\frac{P(A\|B)P(B)}{P(A)}}$	ベイズの定理

加法定理は離散型確率変数での表現になっていますが、連続型確率変数で同時分布から周辺化を行う操作(積分消去)を離散型確率変数に適用したものととらえれば理解しやすいかと思います。あと数式の厳密性の話になりますが、議論の対象の確率変数(連続型 or 離散型)によってPの大文字/小文字を使い分けるようですが、一旦はすべて大文字で記載しました。

最後に書籍を読むために直接必要となった知識ではありませんが、色々と調べていった中でパターン認識と機械学習の応用に使えそうなものをメモしておきます。

名称	URL
線形予測分析	線形予測法 - Wikipedia
ラドン変換	トモグラフィー - Wikipedia

感想と今後

ブログの記事を書くときに間違えていないかなと確認しながら書いているのですが、wikipediaのFFTや畳込みのページがすごく進化していて驚きました。アニメーションが追加されていて直観的に理解しやすくなっています。

Latex初めて使って楽しくなってきたので色々書きましたが、面白い反面なかなか手ごわいですね。 displaystyleつけないとΣの下に式が来なかったりとか。あとは、はてなブログの問題のようですが、行列の&が記述できないので諦めました。

今回はノートの隅にメモしておいただけのものを確認しながら一通り整理して記事にするだけでも、自分の中で色々と整理できました。そろそろ何かプログラム書きながら試したりとかしてみたいですね。では今回はこれくらいで。

2016-05-28

機械学習における回帰と分類について考える

雑記

暑かったり寒かったりと体調を崩しやすい時期ですね。かたりぃなです。機械学習について勉強していますが、まずは基礎的な問題として分類と回帰について考えてみます。 wikipediaや以前買った書籍(PRML)を見ても理論の大半は数式で表されていましたが、ここではできるだけ言葉で整理します。自分の言葉で表現することによって数式の厳密性は失われますが、私自身の理解を深めることになると思っています。特に自分自身が経験してきた仕事の内容で比喩することで考えを整理していくきっかけになると思います。 (そもそも本の内容をただコピペしてしまうと著作権的にも問題になりますし。)

回帰と分類

教師あり学習 - Wikipedia

教師あり学習における回帰と分類はそれぞれ以下のように解釈できます。

入力値と出力値のペアが与えられたとき、この入出力の間にある変換を求めることです。出力値が連続な場合は回帰、離散値の場合は分類になります。そしてこの変換を求めることが学習になります。

とは言っても抽象的でつかみづらいので実際の例をもとにもう少し概念を掘り下げてみます。

回帰

組み込みソフトウェアの開発現場での話ですが、あの分野では性能や品質が重視されます。その背景には

ソフトウェアのアップデートが簡単にはできない
性能要件を満たさないと製品として世に出すことができない
不具合が起きた分野によってはリコールや訴訟、そしてブランドイメージの低下などの大きな被害となりうる

などがあります。インターネットに接続されている機器であればネット経由でのバージョンアップによる修正という余地はありますが、バージョンアップ対象のモジュールによってはアップデート自体がリスクとなりうるケースもあります。具体的な例として、カーネル本体やファイルシステム、周辺機器のデバイスドライバなどをアップデートするには相応のリスクが伴います。 (アップデート中に電源を切られると、ユーザーランド上のアプリケーションならやり直せますが、カーネルやファイルシステムなんかをアップデート中にそんなことされてしまうと、何が起きるかわからない)

私が回帰による分析というものに初めて触れたのがこの品質周りの課題について調べている時でした。ソフトウェアのバグを取り除いていったとしても、一般にハードウェアは経年劣化というものがついてまわります。物理的な動作が含まれるハードウェアはもちろんですし、外見上は物理的な動作を含まないように見えるデバイスも劣化が進みます。フラッシュメモリなどの書き込み回数の上限などがその一例ともいえます。

このような劣化の要因・品質特性を調べるとき、劣化による影響がどのように現れるのかを考えると

劣化の原因は何か
劣化すると何がおきるのか
劣化の速度はどれくらいなのか
工場出荷時の品質はどのくらいのものなのか

などがすぐ思いつきます。これを分析する最も基本的なものが線形回帰による分析です。回帰直線などを求めるアルゴリズムとしては最小二乗法などが有名です。

話を簡単にするために、上記の品質の例をY=aX+b　という中学数学で習う比例のグラフにあてはめられるという前提のもとで考えることにします。劣化の要因は経年劣化であるとします。すると工場出荷時の品質が切片bであり、劣化の速度が傾きaとして表現できます。また劣化の要因は経年劣化なので機器の実稼働時間とします。すると、

ある期間稼働した時点での品質(Y)＝劣化の速度(a) x 実稼働時間(X) + 出荷時の品質(b)

となります。あとはどのくらい劣化したら品質を保てなくなるのか(摩耗劣化であれば、機器から異音がし始める時期など)を推測できます。これで品質特性がわかった！とはできないのが世の中というものです。実際にこれをそのまま使ってしまった場合の問題点として

関係(この例ではY=aX+b)の前提が無いことが多い(測定して分析するまではどのような式が近似に適しているかわからない)
影響している要因(品質工学でいう因子、今回の例では経年劣化)が何であるのかが不明であること

があります。それぞれの問題点について個別に考えていきます。まずは式そのものについて考えると、例えばY=aX²+bという特性かもしれないし、もっと複雑な次数で表される特性かもしれません。品質特性についてはバスタブ曲線と呼ばれるグラフ形状が有名です。

影響している要因について、実際には稼働時間ではなくて電源ON/OFF回数が影響していた場合などが考えられます。要因(品質工学的に因子と呼ばれる)特定するには、デバイスのデータシートにそれら示されていることもありますが、データシート記載の実験環境と実運用環境は往々にして異なるので、そういったものに関してはリスクとコストを天秤にかけたうえで実験・測定する必要があります。しかしながら、劣化に関する実験と測定を行うには多大な時間がかかります。これは摩耗劣化の例でわかるかと思います。

こういった要因を考えるのは非常に難しいもので、実際に測定してみないと何が影響しているか判断しづらい場面というものが多々あります。システムが複雑化している最近のデジタル家電なんか特にそうです。こういった難しさというものはタグチメソッドの田口玄一氏の有名な言葉「技術者の最大欠陥は何が起きるのかわからないと手を打てないということだ」でよく表されるものだと思います。

話が発散しそうなのでここで回帰とは何かをもう一度整理すると、

関係Y=aX+bなどを求めること(そもそもこの式自体を求めることも含めて)が回帰である
求めた式を使って予測することが回帰からの推論である(何年後くらいに壊れることが予想されるか)

と理解できます。製品の品質という課題については品質工学(タグチメソッドや実験計画法など)の分野になるのでここで終わります。ここまでの考えをもとに機械学習の観点から見てみると、

上記の式の右辺にあたるもの(関数)を求めることが学習である(Y=3X+2であるなど)
求めた関数を使って予測を行うことが推論である(Y=3X+2という学習結果があれば、未知の値Xが与えられても目的変数Yを求められる)

といえます。機械学習の分野で学習(学習段階)と決定(推論段階)が分けて考えられている理由も納得できました。

分類

回帰における目的変数Yが離散的な場合が分類と呼ばれる問題です。といっても抽象的すぎるので、具体例で考えます。

身近な例ではゲームにおける当たり判定が分類問題として扱うことができます。話を簡単にするために2Dゲームの矩形で考えます。当たり判定の基本となるアルゴリズムは矩形と矩形の重なった領域があるかどうかです。重なった領域があればヒット、なければノーヒットです。ここで目的変数Yが「当たった、当たっていない」であり、離散値は２値でヒット、ノーヒット(プログラム的表現するならばbool値)です。

単純なゲームであればすべての矩形の当たり判定をチェックすることで目的は十分に達成できます。しかしCAVE(怒首領蜂など)や東方などに代表される弾幕シューティングのように、判定を行う対象物が数百とか数千の敵弾といったことになると少し工夫する必要が出てきます。画面上にある数千の弾丸すべてを毎映像フレームごとにチェックするのは厳しいでしょう。スマホやタブレットのCPUなどでは尚のことです。私がこの問題について考えていた当初はもっと貧弱なCPU環境でした(Pentium-133MHz、SSEはおろかMMXもない時代)。

限られたCPUやメモリ資源の中でこの問題を解決するには、全部の当たり判定を総当たりで処理していくのではなく「当たっているかもしれない」対象物を絞り込みができればよさそうです(そもそも計算する必要すらなくなる)。これはパターン認識や機械学習の分野からみると、ブースティングの一例として考えられます。ブースティングの基本的な考え方は「弱い識別機をたくさんくっつけて(カスケードして)強い識別機をつくれないだろうか」といったもので、顔認識で有名なhaar-like分類器などがこれに該当します。

話を戻して、弾幕の当たり判定を弱い識別機で効率よくやるにはどうするか考えてみます。具体例として、画面サイズは1600x1200, 弾丸の判定は2x2ドット, 自機の食らい判定は4x4ドット、画面上の敵弾の数は10000の一様分布としてみます。総当たりで計算すると10000回の矩形判定をすることになって大変なので、まずは弱い識別機でふるいにかけます。やり方は色々ありますが、ここでは画面全体から分割して考えることにします。まず画面全体に対して自機の位置が右下隅に追い詰められているとしてp(1500,1100)にいるとしたとき、まず左半分にある弾丸には当たるわけがないので、敵弾のX座標値が800以下のものは処理しないことにします。これで判定対象は1/2にまで減りました。 Y座標についても同様にして、判定対象はさらに1/2になります。こうすることで、画面全体の10000個だったものを1/4、つまり2500個にまで絞り込むことができました。これを必要な限り繰り返すことで、当たり判定を行う候補を絞り込んでいきます。ポイントは、最初に大きなふるいにかけるので、自機から遠い弾丸にほど早い段階で判定から取り除かれることになります。

と言ってもこの例ではさほど高速化はできません。もともとが矩形同士の当たり判定という単純な処理なので、絞り込む段階でやっている判定内容は、詳細な当たり判定チェックと同等になってしまうからです。

機械学習から少しそれますが、せっかく振り返ったので書いてみます。じゃあ、どうやったら高速に処理できるか考えていた当初の記憶を振り返ってみると、分岐予測が外れたときのペナルティが大きすぎるので、できるだけ分岐処理(C/C++でいうif)を行わないのが理想でした。結論としては、画面上をマスで分割(16x16とか)して、どこのマスに入っているかを弾丸の移動段階であらかじめ求めておき、当たり判定の段階では自機の入っているマスの弾丸のみ当たり判定を行うといったものでした。たとえば1600x1200の画面を16x16のマスに区切ると、100x75個のマスで画面上を表現できるので、所属するマスを決める処理自体は分岐予測も除算も使わずに記述できるというわけですね。たとえばこんな風に。(xpos,yposともに16bit幅という前提で書いてます。)

x_index = (xpos & 0xFF00) >> 8;
y_index = (ypos & 0xFF00) >> 8;
table[x_index][y_index] = ごにょごにょ;

マスのサイズを16x16にしたのはここでの演算を高速にするためのもので、マスクとシフトだけで済ませられます。機種依存していいなら上記のコードよりもさらに早くxposとyposの上位8bitをテーブル参照のインデックスとして使うだけで済みます。(x86 アセンブラ的にやるとahかalレジスタ読むだけで済んでいた気がします) あとは自機が所属するマスを求めて、そこのマスに入っている弾丸とだけ当たってるかどうか調べれば済みます。自機がマスを跨いでいると面倒だった記憶があります。(とはいっても、この例では最大4マスチェックすれば済みますが)

このマスに区切って考えるという方法は2Dゲームの例では非常に有効ですが、機械学習という観点からは次元の呪いという問題にハマることになります。機械学習では入力データは2Dゲームの例よりももっと次元数が高くなる(たとえば16x16pixelを256次元のベクトルデータとして扱うなど)ので、この方法をそのまま使うことは現実的でなくなってしまいます。

感想と展望

よくよく考えるとこのエントリ自体が回帰分析の一種として考えられそうです。回帰分析をしていた当時は諸事情もあって限られた範囲内での分析しか行えませんでしたが、今やるとしたらもっと色々試せる気がします。具体的には

統計結果がある分布に従うことがわかったとき、出荷台数をもとに数年後の修理台数やコストを予測する
予測されるコストをもとに品質保証にかけるコストとのトレードオフを算出する

などがありそうです。

話ついでに。PRMLの下巻買ってみました。上巻よりさらに手ごわくなった内容ですが色々と興味深い内容ですので、色々試しつつも整理していきたいと思います。上下巻を比較すると上巻が基礎理論、下巻が実用に向けてのさらに発展した内容＋実例といった内容です。

なんか長くなってしまったので今回はこれくらいで。

2016-05-06

色々と本を読んでみた

桜の花も散ってしまって春から夏へと季節の移り変わりを感じます。かたりぃなです。
タイトルどうしようかと考えましたが、記事の内容が整理しきれない気配なのでここで一旦投稿としました。

図書館を歩いていてふと開いた本で新しい世界が開けました。
こういった不意に得られる情報というものが私はとても好きです。電子書籍にはない良さです。
デジタルに比べてアナログな書籍の良いところは

ふとした瞬間に有意義な情報が得られる偶然性(目的以外の一見ノイズとも思える情報の中から得られるものもある)
人間の脳のアナログな記憶にアナログな本はマッチしている(何度も開いてクセがついている、後半のほうといったアナログな参照ができる)

だと思っています。

少し話がそれましたが、今回はニューラルネットワークの続きを少し調べてみたので簡単に整理してみようと思います。
脳科学については専門外なので解釈が乱暴で厳密性・正確性に欠ける部分が多々あると思いますが、ここでは機械学習(そしてディープラーニング)について理解することが目的なので正確性は犠牲にします。私にとっては脳をコンピュータ上で再現することが目的ではありませんので。
ニューラルネットワークの基本概念はこれで一段落にします。

ニューラルネットワークにおける教師無し学習

特徴抽出細胞(feature-extracting-cell)

色々と文献を漁ってみたところ、動物はどうやってものを認識しているのかという生物学的な話にまでさかのぼるようです。
まず「特徴抽出細胞」とはネコを使った実験で発見されたものらしいです。
その昔、ネコの脳の視覚野と呼ばれる部分に特定の傾きのスリット光に反応するニューロンが発見されたそうです。
またスリット光に反応するニューロンは集団をなしていて、ニューロンが反応するスリット光の傾きは連続的に変化していた。と。
なんか日本語おかしいですね。。。
どうやらこれを機械学習の観点から応用したらしく、ポイントは

学習モデルの構造として単純な入力処理だけを行う視細胞と複雑な処理を行う視覚野の２つのモジュールに分けられる
視細胞と視覚野の間にニューラルネットワークを構成
特徴抽出細胞は後天的に形成されたものである
特徴抽出細胞は教師無し学習によって得られたもの

だと理解しました。

この視細胞と視覚野を結ぶニューラルネットワークのことをトポロジカルマッピングというらしいです。
トポロジカルマッピングの代表的なモデルとして

willshaw-malsburg
kohonen

のモデルがあるそうです。

このままではよくわからないままなので、前回の記事のアルファベット認識の例で考えてみます。
機械学習の基本的な概念と歴史を追ってみた - catalinaの備忘録
このときの例の16x16pixel(256次元)のグレースケール画像を使った場合では、我々が画像を認識するうえで当然のように考えているものを落としてしまっていることになります。
もっと具体的にいうと、ピクセル間の距離や並びという概念が抜け落ちることになります。
ここで情報を再構成するための話が、先の生物学の話に出てきた猫のスリット光に反応する視覚野で得られる情報そのものなのですが、視細胞で得られた情報を視覚野で再構成(256次元から二次元に落とす)するための数理モデルのお話のようです。
wikipediaでいうとこのあたりでしょうか。自己組織化写像。
自己組織化写像 - Wikipedia

トポロジカルマッピングの理論の中によく表れるキーワードとしてシナプスの競合と協調がありました。

競合とは周囲のシナプスを抑制しあうこと
協調とは周囲のシナプスを興奮させ合うこと

といった意味合いになります。

最近話題のディープラーニングで使われているという畳み込みニューラルネットワークでは、畳み込み層で二次元カーネルによって画像の局所特徴量を得ることに成功しているようです。
局所特徴量を得られているということは、先にあげた単純なニューラルネットワークでは失われていたピクセル間の影響をうまく抽出できているといえます。
ちなみにカーネルと言ってもOSのコアの話ではなく、画像フィルタで一般的に使われるカーネルのお話です。
具体的なカーネルの種類としてsobelやgaussian、laplacianといえば思い出す人も多いのではないでしょうか。
画像フィルタリング — opencv v2.1 documentation

マルコフのモデル

マルコフ性とは、未来の状態は現在の状態からのみ決定され、過去の挙動とは無関係であるとするものです。
情報処理技術者試験の午前問題で単純マルコフ過程が頻出問題だったと記憶しています。
例えば、
ある日の天気が雨であり、翌日が晴れの確率は30%, 曇りの確率は50%, 雨の確率は20%のとき、２日後の天気が晴れとなる確率は？
などです。
実際の出題では曇りの日や雨の日の翌日の天気がどうなるかといった確率をまとめた表と、これが単純マルコフ過程であると限定したうえでの出題なので、単純な確率計算で答えは求められます。

少し複雑なものでマルコフモデルにおいて直接観測されない状態を含むモデルのことを特に隠れマルコフモデルと呼ぶそうです。

このマルコフのモデルをもとに行うシミュレーションがマルコフ連鎖モンテカルロ法として知られていて、マルコフのモデルの下でシミュレーションを繰り返すと、値が収束していくという特性があります。
情報処理の試験が役に立たないと言われるのは、こういった話にまで踏み込まないで簡単な計算問題や概論だけで止まってしまっているからだと思います。

さて、このマルコフ連鎖モンテカルロ法ですが、何に応用できるのかと疑問に思っていました。
色々と調べてみたところ、後述のベイズのモデルにおける推定のための分布を求めるために利用されるようです。
応用分野においてこの２つがセットで語られることが多いのはこういった関係があったからなんですね。

ベイズのモデル

ベイズの定理を使った推測のモデルです。
認識すべき入力パターンがどの母集団から発生したと考えるのが尤もらしいかをベイズの定理をもとに識別します。
事後分布は事前分布と尤度の積に比例するというのが基本的な考えのようですが、まだ理解が追い付いていないので詳細はまたの機会に。

感想と今後の展望

思いっきり話がぶっとびますが、その昔の有名なゲームにCLANNADというのがありました。人生ですね。
ヒロインの一人、一ノ瀬ことみ(ひらがなみっつでことみ。呼ぶ時はことみちゃん)のご両親は物理学者で、この世界を表す数式を研究しているという設定でした。
もしこの現実世界を理路整然とした数式で表現できるとすれば、ディープラーニングの先にあるものはことみちゃんの両親が研究していた数式なのかもしれなません。夢のある話ですね。
とは言っても、我々人類が置かれている事象の内側だけしか観測できないのと同じように、ディープラーニングも与えられた事象に対してのみ働くだけになるかもしれません。

そんなくだらないことを考えつつも機械学習についてもっと知りたいと思ったので、本を買って読みはじめました。
足りない情報や基礎理論は図書館を利用することにしました。

一冊目。ハズレでした。少なくとも私にとっては。
一通りの理論を理解し終えた人が机の上にリファレンスとして置いておくにはいいのかもしれませんが、参考書としての最初の一冊としては情報不足でした。
www.amazon.co.jp

二冊目。図書館で借りて読みました。
確率の基本理論から入って同時分布、条件付確率、そしてベイズの定理へと理論的で分かりやすい本でした。
全部は理解しきれていないので、また復習するときにお世話になりそうです。
www.amazon.co.jp

三冊目。これは買っておきました。
上下巻セットにするとお小遣いがほとんど消えてしまうので、まずは上巻を買いました。
まずは一通り概要を読んでみましたが良書だと思います。
一回読んだだけで全部理解するのは困難なので、演習問題に挑戦しつつじっくり読み進めたいと思います。
www.amazon.co.jp

その他。
図書館で借りて読みました。
3Dグラフィックの行列とベクトルで詰まった人におすすめしたい本です。
上記のパターン認識の本を読んでいると行列やベクトルの計算が出てくるのですが、正直ほとんど忘れてしまっているのでちょっと時間をとって復習してみました。
この本はただの行列の計算方法だけでなく、基底や直交(一次従属、独立)の話から始まり、ベクトルを異なる基底へと変換する……いわゆる座標変換まで丁寧に解説されていました。
行列を習った学生の時分、何の使うかわからない面倒な計算が本当に苦痛で、便利になるどころか計算の手間が増えただけじゃないかとすら思ってました。
3Dグラフィックの分野では必須ともいえるものですが、この「計算が面倒」という先入観も相まって、当時最新鋭のPentium100MHzのPCで行列の乗算なんて重い処理を避けて、どうやって高速化するかなんて考えてました。
上京してきたばかりの貧乏学生には3Dグラフィックアクセラレータ(今でいうGPU)なんて高級品、手が届くはずもないので。
そんな過去を振り返りつつ楽しく読めました。
www.amazon.co.jp

今回はこれくらいで。

2016-03-28

機械学習の基本的な概念と歴史を追ってみた

雑記

世間ではコンピュータに囲碁、将棋、チェスなどをやらせて人間を破ったなど思考する機械の話題が賑やかですね。
技術的なキーワードでは機械学習やディープラーニングになるのですが、詳しくいことがわからないので色々と調べてみました。
ディープラーニングの基礎となる機械学習について調べてみたので、自分の理解と考えを整理するためにもブログ記事として書いてみようと思います。なのでこのエントリには見落としや勘違いなど、間違いが含まれていると思います。
なお文献によって用語が異なっていたりするので、このブログでも統一がとれていないかもしれません。

着目点

「人間の仕事が機械に取って代わられる」という記事が週刊誌や新聞に掲載されていたりします。本当にそうなのでしょうか？
こうのネタは定期的に出てきていますが、結局のところ機械に置き換えられることによって人間の仕事の質が変化してきただけだと私は思っています。
とはいえ、最近のAI,ディープラーニングなどの技術の中身を知らないまま歴史からの推論だけで決めつけてしまっていてはちょっと勿体ないですし、エンジニアという視点で自分が置かれている環境に影響を与えうる要素に対して意見をもつなら、何かしらの技術的根拠をもとに結論を出したいところです。
今回調べた内容だけでは最近もてはやされているディープラーニングまでたどり着けなかったので、一通り調べてみてそのうえで自分なりの答えを出したいと思います。

そもそも学習とは

まず人間についての学習とは何かという点について、調べてみました。
人間工学などの書籍をあたってみると大きな分類として

自律的学習
帰納的学習
演繹的学習
知識学習

などがあります。

これら学習の概念を整理すると「学習とは"問題解決のために用いる自己の能力"を向上させる手段」
といえます。

ここから安直に考えると「機械学習とは上記の学習を機械で行うこと」と一言になります。しかしそれだけでは私自身の何の学習にもなっていません。
機械学習の歴史をたどって、機械(コンピュータ)に学習させるという点を掘り下げて考えてみます。

機械学習の歴史

図書館で本を借りて色々と調べてみました。
いくつかの文献をあたってみたところ、必ず出てくる有名なものがあるのでここに並べてみます。

神経回路網的アプローチ
1. (ローゼンブラット)パーセプトロン、ニューラルネットワーク
2. (ルーメルハート)バックプロパゲーション
記号論的アプローチ
1. (ウィンストン)アーチの概念の学習
2. (ミッチェル)バージョン空間法
3. (ザデー)ファジィ理論
4. (ハント)Concept Learning system 決定木の学習

それぞれの概要を理解した範囲で簡単に整理してみます

神経回路網アプローチと記号論アプローチ

神経回路網アプローチとは、最近流行のディープラーニングもそうですが、神経網モデルをコンピュータ上で表現し、その神経網に学習をさせようというアプローチです。
対する記号論アプローチは、述語表現やルールとして記述されたモデルによって学習をさせようというアプローチです。
またこれとは別にPAC学習というアプローチもあるようですが、今回の資料だけでは詳細情報は得られませんでした。

神経回路網アプローチ

パーセプトロンとニューラルネットワークそしてバックプロパゲーション

パーセプトロンは一つの神経細胞を表します。パーセプトロンは複数の入力を受け、複数の出力を行います。
パーセプトロンの入力は入力値に対して重みを掛け、入力の合計値にバイアス値が足されます。これを活性化関数(古くからの実装ではシグモイド関数)に受け渡すことで出力が決定されます。
ここで学習結果は重みの値となります。
このパーセプトロンを組み合わせたものがニューラルネットワーク(NN)と呼ばれるものです。
NNには2種類あって、階層型と相互結合型(再帰型)に大別されます。
現在広く使われているものは階層型NNのようです。
階層型では層の呼び方としてN層ニューラルネットワークなどと呼ばれます。
このとき、入力層、中間層(隠れ層)、出力層の３つの呼び方で入出力が区別されます。
入力にも出力にもあたらない層が隠れ層です。

さて、パーセプトロンに話を戻して、
先にあげたアルゴリズムのままでは何がすごいのかイマイチわからないので、次のような例を考えます。

2入力1出力のパーセプトロンで、2入力の重みをそれぞれ1とし、バイアスを1.5とした場合
2入力1出力のパーセプトロンで、2入力の重みをそれぞれ1とし、バイアスを0.5とした場合
1入力1出力のパーセプトロンで、入力の重みを-1, バイアスを-0.5とした場合

手で計算して真理値表を書くと見覚えのある表になるのですが、このパーセプトロンはそれぞれ、AND,OR,NOTの論理回路になります。
これは重み値とバイアス値がある値をとるとき、基本論理回路(ゲート)として振る舞うことができるということを意味します。
すなわちコンピュータが学習によってコンピュータの基本素子ともいえる論理回路を構成できてしまうということです。
あとは組み合わせ次第でどんな回路も作れますし、加算器や乗算器、シフト回路、果ては記憶素子としてのFF(フリップフロップ)やラッチを作ることだってできそうです。
機械学習の結果で論理回路が構成されるなら、機械学習の結果次第で何でもできる！！
なんとなく機械学習の神経網アプローチってすごい気がしてきました。

機械学習としてニューラルネットワークを利用するためには、学習の結果をパーセプトロンに反映する必要があります。
このために教師アリ学習として扱ったとき、「NNから得られた出力信号」と「期待する結果」とを比較することで誤差を知ることができます。
この誤差をもとに重みを調整すればよいのですが、この計算方法はローゼンブラット式とヘブの学習と呼ばれる式で違っていました。
学習結果の反映のためのパラメータ調整の話になってくると思うので、ここでは複数の方式があるということだけ記しておくことに留めます。

バックプロパゲーションですが、乱暴に概要だけを述べるならば、「学習結果を出力パーセプトロンから入力方向へどんどん遡っていき、中間層(入力でも出力でもない層。隠れ層とも呼ばれる)の誤差を小さくしていく(学習を行う)こと。」といえます。
このバックプロパゲーションがどうやら革命的なことらしく、この技術が登場するまではニューラルネットワークの研究の冬の時代と呼ばれていたらしい。
この技術が登場したおかげで、パーセプトロンを多層にできるようになったとのことです。

NNの実際の応用事例ですが、たとえば画像からのパターン認識などに利用できるようです。
簡単な例では16x16ピクセルのグレースケール画像に示されるアルファベットの認識。
この例では入力は16x16=256次元のベクトルとし、入力層に割り付けます(入力層は256個のパーセプトロン)。
出力層は26種類のアルファベットのいずれか１つなので出力層の26個のパーセプトロンのうち最もよく反応したものを識別結果とする。
などがあります。

私の中でもまだ理解が追い付いていない＆疑問が残ったままの要素として、次の２つがあります。

NNでは過学習が問題となったようだが、これに対するフォローはどうするのだろう
隠れ層の数はどのようにして決めるのが適切なのだろう

まだ概論しか触っていないので、これらの問題はもっと深い部分で解決しているのかもしれません。

ニューラルネットワークに対する個人的見解

NNは今後もどんどんと進化する可能性はありそうです。
相互結合型(再帰型)のほうも興味深いのですが、今回あたってみた資料では情報が古いせいか相互結合型NNの情報が少ない状態でした。相互結合型NNは結果の収束の判別方法などまだ理解しきれていない点が多いですが、個人的にはすごく興味深い話だと思います。
NNをゼロから実装するのは手間がかかりすぎるので簡単に実験できないかと調べてみたところ、例えばopencvのmlモジュールに含まれていますし、R言語でもNNetパッケージが提供されています。
opencvのmlモジュールでは多層NNが実現できますが、学習結果を可視化する手間がちょっとかかりそうなので、試してはいません。
R言語のNNetパッケージは3層(隠れ層が1個のみ)のNNしか使えないのがちょっと残念ですが、学習結果のNNを可視化するパッケージが某大学の某教授が公開してくれていますので、動作概念を確認するのに利用できました。
R言語での多層NNパッケージも調べればあるのかもしれませんが、それはまたの機会ということで。
NNを試してみて思ったことは、動作検証の方法が統計的手法になりがちなので、いままでやってこなかった統計的手法についても調べていく必要がありそうです。

記号論的アプローチ

ウィンストンのアーチ

「積み木のパーツの中から３つだけ使ってアーチを作れと」言われれば、ブログ書いたりネット見たりしている我々ならすぐにできるものですね。
技術の無駄遣いをする人ならもっととんでもないものをつくってくれるかもしれませんが。
アーチの回答ですが、人間がやるとすれば次のような解が導き出されます

柱を2個立てる
柱の上に何かを1個載せる

これを使って「コンピュータにアーチの概念を学習させる」という論文が「ウィンストンのアーチの概念の学習」と呼ばれるものです。

概念は非常に分かりやすく直観的なものでした。
直方体や円柱、三角などの積み木を意味ネットワークと呼ばれるグラフで表現し、ノードを積み木の属性、ノード間のつながりで積み木の関係を表現します。
このグラフがアーチの概念を表すまで正の例、負の例を与えて学習をさせるというものです。
興味深いのは、神経網的アプローチで神経細胞の出力を数値で表すのに対し、こちらはmust,notなどの論理式で表すということでした。
負の例の中でもニアミスは概念を学習するうえで非常に重要な位置をしめているとのこと。
上記の人間がやる例でも見落としになっていますが「柱同士がくっついていてはダメ」という概念を学習によって得られるとしています。
柱同士がくっついている負の例(touchというノード間のつながりをデータとして与える)を与え、帰納的に学習させるというものでした。

ウィンストンのアーチ自体は帰納的にコンピュータに概念を学習させるという良い例だと思います。
この方法論をそのまま再現しても意味はありませんが、論理的な学習手順のアプローチの基本原理として非常に参考になる例でした。
問題点もわかっていて、ニアミスは良い学習になるが、学習済みモデルとの差が1要素のみであることなど制約は大きいです。
この理論を発展させたものが決定木であり、決定木を弱識別機として多数組み合わせたものがランダムフォレストです。

決定木

決定木は与えられた訓練事例から、事例の正負を判断するための木を使ったアルゴリズムです。
アルゴリズムとしては

ハントの Concept Learning system
クィンランの ID3
C4.5

こちらのアルゴリズムも帰納的に概念を学習させるモデルで、画像処理ライブラリなどに含まれていたりします。

決定木は与えられた事例集から、それを識別するための概念による分類を行うために利用できます。
例えば、ある有名キャラクターを概念的に識別することを考えると

性別(o女性,x男性)
髪の色(o金,xピンク,x緑)
髪の長さ(oロング,xショート,xなし)
肌の色(o肌色,x茶褐色)

などと識別するための木を作ることができます。これが概念の木(=決定木)そのものになります。
先のキャラクターの例で考えたとき、あるキャラクターそのものを言い当てられるのは当然として、そのキャラクターと同様の特徴をもつ母集団全体の性質を決定木は示していることになります(女性、金髪、ロング、肌色)。
また、決定木の高速化に大きく寄与していると考えられるアルゴリズムとして「エントロピーや情報獲得率を考慮した木の構築」が挙げられます。
これは特徴を言い当てて答えを絞り込むために、「何を優先的にチェックすべきか」を学習によって獲得していることに相当します。
エントロピーとその正規化について考えると、学習のデータセットにデータベースでいう主キーに相当するものが含まれていたりすると、その主キーのみで全てが決定できるため、情報量が最大となってしまい、期待した決定木が生成できません。
言葉でいうと「場合の数」が多い概念(髪の色を500種類に分類したなど)が優先されてしまうということですね。
(これは学習データセットに対する識別率は100%だが、汎用性は皆無な木が生成されてしまうことになります。)
そのため、エントロピーの正規化した情報獲得率とも呼ばれる方法が利用されるようです。

ファジィ理論

バブル崩壊と叫ばれていたころ、家電製品でどこもかしこもファジィファジィと謳っていたのを記憶しています。
父に聞くと「ふつうの機械はON/OFFしかできないけど、ファジィなら「いい塩梅」に調整してくれるんだよ」と答えてくれたのを覚えています。
今回図書館で借りた本だけでは詳しい情報は得られませんでしたが、父の教えてくれた答えは「だいたい合っている」気がします。

ファジィ理論ではこのような「だいたい合っている」が示す範囲を「ファジィ集合」と「メンバシップ関数」で特性づけた集合で表すそうです。
メンバシップ関数がその集合に属す度合を示します。
ファジィ集合は集合演算が定義されていますが、一般的に呼ばれる集合(クリスプ集合)とは若干異なるものです。

今回あたった資料では概念だけしか紹介されていませんでしたので、概念イメージのみに留めます。
私が理解した範囲でのファジィ制御の概念として、
あるアナログ入力値に対して上記の「ON/OFFしかできない」と「いい塩梅」の制御を行ったときの概念イメージを表現したとき

ON/OFFのみでは、制御範囲は理想的なステップ応答のような矩形になる
いい塩梅では、目標値を中央とした分布(メンバシップ関数で変わる)に属するファジィ集合に応じた制御

となりそうです。
実際の家電制御では計算量が増えてしまう分布からの計算ではなく三角形近似などで十分であるケースがほとんどのようです。

ファジィ理論とニューラルネットワーク

人工知能の概論と歴史を駆け足で追ってみましたが、ファジィ理論とニューラルネットワークを組み合わせることでもっと効率的な学習が行えるかもしれませんね。
と思ってネットで調べてみたらありました。ファジィ・ニューラルネットワークというらしいです。まんまですね。

感想・今後の展望

今まで理屈を知らなかったものを知るということは非常に楽しいものです。
また、新しい知識を得てもまたわからないことが出てくるのでそれはそれで楽しいものです。
最近話題の機械学習はもっぱらニューラルネットワークが基礎的技術となっていますが、試しにR言語でnnetパッケージ動かしてみたところ、結果がわかりにくすぎる状態でした。(学習結果のネットワークを見てもそれが期待する結果なのかどうか判断できない)
学習結果の尤もらしさを人間が客観的かつ迅速に判断する仕組みづくりがどこかで必要なんじゃないかなぁと思います。
それがいわゆる統計的アプローチにつながっていくのだろうと予想しています。
神経回路網的アプローチに対する記号論的アプローチでは述語表現や概念グラフが学習結果として現れるので、学習結果を人間が確認するのも簡単そうです。
しかしながら、今回は単純な例しか考えませんでしたが、複雑な学習になってくると学習結果の確認すること自体が現実的ではない巨大な概念グラフになってしまうでしょう。
色々と調べて考えてみた結果、ビッグデータやデータマイニングと騒がれたときと類似した問題解決手法とそれぞれの課題が見えてきました。
また、根底にある考え方として2種類(記号アプローチと神経アプローチ)のアプローチがありましたが、それぞれアプローチは異なるものの、どこかで一本の糸につながりそうでつながらない状態です。
まだ考えの整理が追い付いていない部分が多々ありますが、自分なりに考えてみて、色々やってみたいと思います。

今回はアルゴリズムや機械学習や人工知能の仕組みと歴史を中心に調べましたが、次回は少し趣向を変えて統計的アプローチ周りを見てみたいと思います。
具体的には統計的決定論と呼ばれるものですね。ベイズとかマルコフとか。

2016-03-06

Rでコスプレ用スカートの可視化

R 雑記

結城先生の数学ガールを読み始めてみました。かたりぃなす。
もっとゆるい系の本かと思っていましたが、中身はすごく真面目かつ分かりやすく説明されていて、良書だと感じました。
20年近く前にフーリエの冒険を読んで以来のよい数学書です。
こういう本は非常に良い刺激になります。

そんなわけで、少しばかり数学的なことをやってみようと思いました。

コスプレ用のスカートの寸法の自動計算と可視化

コスプレ用のスカートの型紙についてですが、これが結構色々な局面で使います。
魔導士系のローブなども全円スカートを応用して広がりを少なくすることでそれっぽく作れたりもするので。
有名どころのサイトではこちらで簡単な計算ができます。
http://yousai.net/nui/seizu/furea/zenengomu.htm

全円スカートなら丈を設定して、出てきた寸法をそのまま使えば事足ります。
ただ、先にあげたように応用して作るときは型紙の書き直しが必要だったりするので、出来上がりがイメージしづらいです。
ロング丈のローブやスカートになると型紙作るだけでも大変なので実験するにも躊躇してしまいます。

よくよく考えると、衣装を作るときは毎回紙に書いて似たような計算しています。
ここはエンジニアらしく解決してみたいところです。
数学ガールも読んで気分も上がってきています。
そうだ、久しぶりにR言語を叩こう。
データの可視化だ。
というわけで作ってみました。

私の型紙づくりのやり方で特によく使う数学的なものとして、

三角形の合同・相似
円錐・円柱・球などの空間図形
一次関数に関する基礎的なもの
円の公式に関するもの
垂直二等分線、平行線などを利用した作図

などです。
義務教育で習う内容がほとんどですね。

実際にどんな局面で利用するのかというと、すぐ思い浮かぶものとして

スカートの丈を調整するとき
- 円錐として考える
  - 空間図形
  - 円の公式一般
  - 直角三角形の一般的なもの
製図するとき
- 垂直二等分線を使って正確な垂直線を書く
- 大きな円弧を描くとき、三角形からの近似(歪まずに書ける程度まで三角形分割する)

などです。

今回作ったプログラムは次の3つのパラメータをもとにどんなスカートができるのか可視化します。

スカート丈(パニエなどで広げている前提)
ヒップサイズ
着用時のスカート裾部分の半径(広がり具合)

ヒップサイズ80cm, 丈100cmのロングスカートを可視化する実験です。
見栄えに一番大きく影響するパラメータはスカート半径なので、このパラメータを変動させてどうなるか見てみることにします。

スカート半径40cm

それほど広がらない魔導士向けローブなどで使えそうです。
女性向けなら丈を短くしてタイトスカートなどに使えるかもしれません。
f:id:Catalina1344:20160306154322p:plain

スカート半径70cm

ワイヤーパニエなどでスカートを広げたロリータ系の衣装やドレスでよく見かける形状です。
これはいわゆるAラインと呼ばれているものですが、個人的にはプリンセスラインのほうが好みです。
R言語なら数式回りの機能は豊富なので、工夫すればプリンセスラインの曲線くらい簡単に表現できる気がします。
しかし疲れたので今日はこれくらいにしておきます。
f:id:Catalina1344:20160306154331p:plain

読みづらいですがコードはこちら。
とりあえず書いてみたレベルなので間違いあるかもしれません。

# ヒップサイズ(cm単位)
west_size <- 80
# パニエ着用のうえ最大限に広がったスカートの半径(cm単位)
skirt_diameter <- 50
# スカートの丈(cm単位)
skirt_length <-	100

# ヒップサイズから円錐上部の円の半径を求める
top_R <- west_size / (2*pi)

# 三平方の定理を使ってスカートの高さを求める
# スカート丈=斜辺, X=スカート半径-ヒップサイズ, Y=hight
dx <- skirt_diameter - top_R
hight <- sqrt(skirt_length^2 - dx^2)	# 実際に生地を使用するY軸の寸法

# スカートの半径(xy平面)とスカートのY寸法から、中心からの距離に対するZの増分を求める
deltaZ <- hight / dx

# Zの増分をもとに、円錐の頂点を求める(x=0,y=0のときのZの値)
y_max <- deltaZ * skirt_diameter

# ヒップラインを求める
# 円錐の頂上はヒップラインでカットしてスカートらしく見せたいので
hip_line <- y_max - hight

# 以下はグラフを書くための処理
# 3Dグラフ用のX,Y軸のデータを準備
x <- seq(-skirt_diameter, skirt_diameter)
y <- x

# z値を求めて3Dグラフ化用のデータセットを準備
zh_func <- function(x,y){
	xy_2dvec <- sqrt(x^2 + y^2)	# XY平面における中心からの距離
	xy_2dvec * -deltaZ		# 中心から遠ざかるほどスカートは垂れ下がる
}
z <- outer(x,y,zh_func)

# 最小値,最大値を超える値の切り捨て(ヒップラインとスカートの裾)
threshold <- function(v){
	if(v > -hip_line)	return (-hip_line)
	if(v < -y_max)		return (-y_max)
	else			return (v)
}
z2 <- apply(z, c(1,2), threshold)


png("output.png")

# 三次元グラフ出力
v_theta <- 10
v_phi <- 40
persp(x,y,z2, 
	scale=FALSE,
	theta=v_theta, phi=v_phi,ticktype="detail", shade=0.5)

dev.off()

感想

せっかくのR言語なのに統計的なことは何もやっていないですね。
今回はただ可視化してみただけなので、三次元グラフのスケールが統一できていないのがちょっと残念です
グラフのスケールが毎回異なってしまうと、たとえばスカートの丈だけが異なるようなパラメータでグラフを生成したとき、見た目からは差が分かりにくいです(軸の値を読まなければいけない)。
まずは可視化できたのは一つの成果なので、良しとします。
色々と思うところはありますが、今は自分が楽しいと思えることに打ち込みたいですね。

今後の展望

これをそのまま型紙に起こせるようになれば色々と便利に使えそうな気がします。
今回のような全円スカートであれば、どこか一か所をカットしてそこから展開できるかもしれません。
プロットした座標間の距離を保ったまま平面上に落とし込んであげればそれっぽく展開できるだろうとは思うのですが…。
また気が向いたら試してみます。
あとコスプレ用衣装にありがちなフリルやギャザーを入れた場合も同じように可視化できると楽しそうです。

2016-03-03

ARの原理実験

OpenCV openGL

ARの原理について色々と調べてみました。
結論からいうと、ライブラリを２個組み合わせるだけでそれらしいことはできそうです。
この原理試験を簡単ではありますが試してみたいと思います。
世間ではOpenGVとかvuforiaとかAR向けのフレームワークが出てきていますが、
今回は原理試験が目的なので、フレームワークは使わずにOpenCVとOpenGLを使ってみたいと思います。

ARを実現するための基礎技術として

カメラ内部パラメータの推定
カメラ外部パラメータの推定(物体の姿勢推定)

が大きな課題だと思います。
検出対象物の姿勢と位置さえわかれば、あとは3Dグラフィックエンジンを叩いてレンダリングするだけです。
今回はこの２つの基礎技術について実験をしてみました。

実験の目的

今回の実験の目的ですが、
カメラキャリブレーション(カメラ内部パラメータ推定)を行わずに手軽にARを実現できないか？
と思ったのが事の発端です。
例えばARアプリを作ったとしてもキャリブレーション用の器具(チェスボードとか)をユーザーに持ち歩けというのはあまりに無茶ですし、
何かしらの既知の大きさをもった物体(例えば後述のカード)を机に並べるとかでできなくはないですが、手間がかかってしまいます。
色々と調べていくうちに、そもそもキャリブレーションしなくてもいいんじゃないか？という結論に至ったので、
こんな理屈で本当にARらしいことができるのか？を確認することが目的です。

理論的なもの

ARにおける物体の姿勢推定の前に基本的なカメラモデルの理論です。
多くの場合では問題を簡単にするためにピンホールカメラモデルが利用されているようです。
ピンホールカメラモデルは次のような式で表されます。
式：撮影画像上の点 = 内部パラメータ行列 x 外部パラメータ行列 x マーカー座標系の点
論文やwebサイトをあちこち見て回りましたが、最終的にはこのような式で表現されます。
この式には含まれていませんが、実際には歪み補正パラメータというものもあるのですが、次の2つの理由から無視することにしました

広角レンズでも使っていない限りはそれほど大きくない
最近の一般的なwebカメラ(特に手元で実験に使うwebカメラ)はそこまで歪みがない

なお、OpenCVで提供されているカメラモデルも同様のものです。
カメラキャリブレーションと3次元再構成 — opencv 2.2 documentation

上記式のパラメータの詳細は論文をあたってもらったほうが正確な情報が得られるので、ここでは深く説明はしません。
(私自身も突っ込んだ質問されても正確に答えられる自信がないというのもあります。)
この式で示されている各パラメータをソフトウェア上で求める方法について軽く見ていきます。
「内部パラメータ行列」ですが、これを求める(推定する)行為がカメラキャリブレーションと呼ばれています。
「撮影画像上の点」はカメラから取り込んだ映像フレームに対して物体検出などの画像処理を行うことで求めることができます。
「マーカー座標系の点」も同様です。

残ったの「外部パラメータ行列」は今回のエントリの主な着目点です。
マーカー座標系の点と撮影画像上の点を正しく対応付けられているとき、上記式のうち未知である外部パラメータ行列を推定することができます。
この行為をARの世界ではマーカーの姿勢推定などと呼ばれています。(アルゴリズム的にはPNP問題を解くとも呼ばれます)
ここでのPNP問題とは「Perspective-n-Pont Problem」のことで、P≠NP問題のことではないです。

今回は平面マーカーを利用するので、全ての点が平面上にあるという前提のもとに姿勢推定を行います。
そのためには4つの対応点が必要とされるようですが、まだ詳しい理屈がわかっていないです。
この点についても追々調べていきたいですね。

さて、画像処理に対するお話は一旦おいて、3Dグラフィックスの概念レベルのお話です。
3Dグラフィックスでジオメトリをレンダリングするための頂点情報の変換手順は概ね
クリッピング座標系 = 射影変換行列 x ビュートランスフォーム x ワールドトランスフォーム x ジオメトリ座標の点
となっています。
実際のレンダリングパイプラインではこのクリッピング座標系から正規化デバイス座標系を経て、スクリーン座標へと変換されます。

この3Dグラフィックレンダリングの手順とカメラモデル、どちらも「スクリーン投影までの変換手順」を定義しているといえます。

座標変換の手順は表現方法や途中で経由する座標系こそ違いますが、
ピンホールカメラモデルも3Dレンダリングパイプラインも目的は同じです。

本題のARですが、ARとは「実世界から取り込んだ映像上に任意の3Dモデルをオーバーレイする」と考えたとき、カメラモデルとレンダリングパイプラインで同一の変換が行えるようになっていれば、話は簡単になりそうです。
具体的には
映像入力について

実世界の立体像->変換(ピンホールカメラモデル)->ソフトウェア内で表現された画像

映像出力について

ソフトウェア内の3Dモデル->変換(3Dレンダリング)->オーバーレイ表示する画像

の「変換」にあたる部分が同一の変換を行うのであれば、簡単にARが実現できそうです。

数式的な証明はできていませんが、おいおいやっていきたいと思います。(ARには魅力を感じているので、基礎は押さえておきたい。)
大昔に3Dグラフィックスやり始めた頃にスクリーン投影に関する部分は深くは追っていなかったのが悔やまれます。
やっぱり基礎って大事ですね。

さて、ここまでの考えが正しいとした場合、

3Dグラフィックスエンジンを使ってキャリブレーション器具をレンダリングし、そのレンダリング結果からカメラパラメータを推定できるのではないか？
3Dレンダリング結果をもってしてキャリブレーション可能であるならば、実カメラを使ったキャリブレーションは不要ではないか？
さらに突き詰めると、レンダリングパラメータが一意に決まるのであれば、キャリブレーションの手順すら不要で、内部パラメータも一意に計算できるのではないか

などが考えられます。

実世界のカメラをプレビューのみに利用することになるので多少の誤差が出ることは想定されますが、エンターテインメントとしてARアプリを作るなら多少の誤差であれば許容されそうです。(手軽に使えるほうが大事)
逆にミッションクリティカルな分野ではもっと正確なアプローチが求められると思っています。

というわけで本題の実験です。
姿勢推定の実験が目的なので、まずはそれ以外の部分は適当にでっちあげます。

物体検出

本格的に作るのであれば適切な特徴量の検出器などの採用などを検討するところですが、今回の実験でそこまですると本題から逸れてしまいます。
というわけで、実験環境を限定して目的の機能を作ることにします。

カメラは手元のwebカメラ
机の真上から、真下に向かって撮影する
検出対象は机の上に転がってたMTG(MagicTheGatharing)の黒枠カード
ノイズが邪魔なので、白の厚紙の上にカードを置く

この限定条件下であれば、次のような簡易的な手順でカードを識別できます

カードのイメージを事前にスキャナで取り込む
入力映像フレームに対する事前処理
輪郭抽出
スキャナ取り込み画像と入力映像、それぞれの輪郭のモーメントの比較

若干の誤検出はありますが、とりあえずの実験としては目的を充分に果たせます。

物体の点を求める

物体の輪郭のままではマーカー座標系の点との対応付けが困難なので、適当にそれっぽくなるようにしました

輪郭を長方形に近似し、4つの点を求める
4つの点をマーカー座標系と対応付けるため、並び順を決める(右上から時計回りに4つ)

これに対応するマーカー座標上の点は物理的に対処します。

道具箱を開ける
ノギスもしくは物差しなどの測定器を取り出す
測定(物理)
物体の点と同様の並び順で座標をソースコード中に定義する

これでマーカーとの点の対応付けもできました。(画像中の座標、物理座標どちらも右上から時計回りに4点あります。)
とりあえずの実験らしく超乱暴な逃げ方ですね。

やっと実験の目的のところまで来ました。

OpenGLでレンダリングした結果を使ってキャリブレーションする

キャリブレーションの関数自体はOpenCVに用意されているのでそれを使います。
方法は何でもいいのですが、よく使われているZhangの手法を使います。

チェスボードのイメージを作る
チェスボードのイメージをテクスチャとしてOpenGLに転送する
板ポリゴンに上記テクスチャを貼ったものをレンダリングして、その結果を取り出す
取り出したレンダリング結果の画像を使ってキャリブレーションを行う

ここまでで、以下のパラメータがすべてそろいました。

撮影画像上の物体の点(物体検出)
カメラ内部パラメータ(キャリブレーション)
マーカー上の点(物理的に測定)

カメラからの入力画像中に含まれる物体の姿勢推定を行う

これらを使ってマーカーの位置・姿勢推定を行います。
OpenCVの関数でいうとsolvePnpです。歪み補正パラメータは無視したいのでcv::Mat()を与えておけばいいです。
推定結果が正しいかをとりあえず確認してみると、どうやらそれっぽい推定はできています。
こういったものの確認にはOpenCV3.0のarucoにあるDrawaxis関数がとても便利です。

推定した姿勢情報、位置情報をもとに3Dオブジェクトをレンダリングする

あとは簡単！
とか思っていたら足元をすくわれました。
しくじったポイントをメモに残しておきます。

回転の定義がOpenCVとOpenGLで違う

solvePnPが返してくる姿勢推定結果の回転の情報は、xyzの回転ベクトルであり、ベクトルの向きが回転軸、ベクトルの大きさが回転量を表します。
対するopenGLのglRotatefなんかは回転軸と角度を度数で指定します。
変換方法を考えましたが、ちょっと面倒ですね。

そもそも私がやりたいのは姿勢推定で得られた結果をOpenGLでのレンダリングに反映したいだけです。
OpenGL1.x系のAPIをわざわざ叩く必要はないということに気付いたので、glLoadMatrixで行列を直接投げ込むことにしました。
回転行列の生成はopenCVにロドリゲスの公式によるものが実装されているので、コレを使います。
これで3x3の回転行列が生成できました。
行列が完成したら、OpenCVからopenGLの行列表現に変換するために転置してあげればOKです。
ついでに、OpenGLの行列は4x4なので、定義に従って残りの成分に平行移動成分などを設定します。

OpenCVとOpenGLで座標系と行列の定義が違う

OpenCVのカメラキャプチャの映像をそのままOpenGLのレンダリングバッファに転送してしまうと上下反転します。
座標系違いますもんね。。。

OpenCVはX+は画面右方向、Y+は画面下方向
OpenGLはX+は画面右方向、Y+は画面上方向

です。
カメラからの入力画像自体を毎フレーム上下反転するのはコスト高なので、OpenGLのカメラ行列で解決します。
カメラ行列の生成時にカメラの上向き方向を定義するとき、Y+が画面下方向になるように指定します。
カメラを逆さまにするってことですね。

カメラからの入力画像しかレンダリングされない

実験ではOpenGL1.x系を使っているのですが、レンダリングバッファに対してカメラからの映像を書き込んでしまうと、深度バッファも更新されてしまいます。
そんな深度バッファのところに3Dオブジェクトをレンダリングを試みても期待した結果は得られません。
カメラ映像を書き込んだ後で深度バッファのみクリアしておく必要があります。
ARでは奥から順に

現実世界の映像
オーバーレイする3Dオブジェクト

となるので、普通のARの使い方であればこれで事足りそうです。

魔法っぽいものをレンダリングしてみた

3Dモデリングは苦手なので、なんとかARらしいことをしてみました。
ARって魔法みたいですよね。カードゲームの世界も剣と魔法のファンタジーですよね。
それっぽくなるようにレンダリングしてみました。

板ポリゴンに六芒星を描いたα付きイメージ(inkscapeで作った)
パーティクル用のフレア画像α付き(gimpのグラデーションフレアで作った)
六芒星、パーティクルともに加算半透明
目立つように画面全体の輝度レベルを半分に落としています

f:id:Catalina1344:20160303172057p:plain

静止画なので分かりにくいですが、カードの姿勢推定も動いているので、傾けたり回転させたりするとそれに反応します。
(物体検出がただの輪郭抽出アルゴリズムのせいで、カードを直接手で持ってしまうと検出できなくなります。姿勢は台紙を傾けて確認しました。)

感想・今後の展望

妙なところで躓きましたが、今回の手法であればそれっぽくARできました。
まだ懸念点・疑問点が残っているので、ここに記しておきます

使用したwebカメラは固定焦点だが、オートフォーカスなカメラではこの手法は通用するのだろうか？
solvePnp(PNP問題を解く関数)の具体的なアルゴリズムと特性を知っておきたい
レンダリングパラメータから一意にカメラ内部パラメータを算出できるのではないか？
手などで遮蔽されるとカードが検出できなくなるのはかっこ悪いので、なんとかしたい。

キャリブレーション結果のカメラ内部パラメータをメモしておきます。しっかりと理論を理解してから計算で求めたほうがよさそうですね。

[580.5345906650882, 0, 320.354109100925;
 0, 580.5217311914437, 239.9422018792191;
 0, 0, 1]

参考にさせていただいたサイト、文献など

http://qiita.com/akira108/items/e5b43810b64cd921a947
http://www.cyber.t.u-tokyo.ac.jp/~tani/class/mech_enshu/enshu2011mi2.pdf
http://www.wakayama-u.ac.jp/~chen/education/cv/pinpara.pdf

2016-02-01

TWE-Liteで外部センサとI2C通信

電子工作

ずいぶんと長い間放置してたものを再開してみました。かたりぃなです。
再開したのはこれ。
catalina1344.hatenablog.jp
もう一年くらい前なんですね。

構成の再検討

TWE-LiteのライターがないからArduinoでなんとかできないかなーとかやっていましたが、
電圧が5Vと3.3Vで違っていたり、色々面倒になってきたのでTWE-Liteのライター買ってきました。
これでTWE-Liteとセンサ、LEDだけで目的達成できそうです。
http://akizukidenshi.com/catalog/g/gM-08264/

TWE-Lite Rの工夫ポイント

そのままではちょっと使いづらかったので、写真のように二階建てにしてみました。
これでブレッドボード上でなんでもできますし、USB給電の電流を気にする必要もほぼなくなります。
(ジャンパーピンでUSB給電OFFに設定して、外部から電源引っ張ってきてます)
f:id:Catalina1344:20160201220115j:plain

センサのアナログ値を直接TWE-Liteに突っ込んでみる

使うセンサは三軸加速度センサKXP84にしました。(道具箱に転がってた)
http://akizukidenshi.com/download/KXP84-2050_Specifications_Rev2.pdf
このコのアナログ出力のXYZの電圧をそのままTwe-Liteのアナログ入力に突っ込んでみましたが、ダメでした。
KXP84のレンジが1.2V～3.5Vであるのに対して、Twe-LiteのADCは0～2.4Vのようです。
アナログの世界では太刀打ちできないのでデジタル世界でなんとかできないかと考えてみました。
TWE-LiteではI2C使えるのでこいつでなんとかなりそうです。
センサ側はSPIもついてるのですが、TWE-Lite側に実装がなさそうなので、I2Cメインで考えます。

TWE-Liteのコードの全体を眺めてみる

東京コスモス様から公開されているソースコードを落としてきて全体を眺めてみました。
やってることは組み込みソフトウェア的だなぁといったところです。(当たり前か)
組み込み系特有な部分はありますが、今回の目的だけであればふつうにC言語でペリフェラルのライブラリ叩いてるだけなので、そこまで難しくはないと思います。
マニュアルも公開されているので、読みながらやれば目的の機能を追加するくらいはできそうです。(その結果、全体の整合性がどうなるかは未確認)
「超簡単！TWEアプリ」はやってることが多いので読むなら機能を絞っておかないと迷いそうです。
特に低消費電力向けにスリープモードがいくつもあったりとか。このあたりはたぶんZigBee固有なお話だったはず。
http://mono-wireless.com/jp/products/TWE-ZERO/index.html

I2Cの実装を眺めてみる

BH1715というチップとI2C通信する機能があったので眺めてみました。
センサの値が確定するまでタイマイベントを待つ以外はふつうのI2Cの通信のように見えます。
I2Cの仕様詳しくは知らなかったのですが、BH1715の資料が日本語だったこともあり、コードと照らし合わせやすかったです。
http://rohmfs.rohm.com/jp/products/databook/datasheet/ic/sensor/light/bh1715fvc-j.pdf?_ga=1.179572772.1693979872.1454157489

マネして書いてみたコードがこちら。
I2CプロジェクトのMain.cにあるvHandleSerialInputで適当なキー入力をトリガとして以下のコードを起動するようにしました。
BH1715と違ってセンサの測定結果が出るまで待たなくていいので単純です。

    bool_t bres = bKXP84startRead();
    if (bres) {
        vfPrintf(&sSerStream, LB "KXP 84 write success");
    } else {
        vfPrintf(&sSerStream, LB "KXP 84 is not found.");
    }
    // Sr + SAD + R + RA
    vfPrintf(&sSerStream, "I2C Bus read start");
    uint8* recv_data = i16KXP84readResult();
    vfPrintf(&sSerStream, "KXP84.read finish");
    vfPrintf(&sSerStream, "X=[%X%X],Y=[%X%X],Z=[%X%X]",
        recv_data[0],
        recv_data[1],
        recv_data[2],
        recv_data[3],
        recv_data[4],
        recv_data[5]);

ここで呼び出しているbKXP84startReadとi16KXP84readResultの関数もマネしながら書いてみました。

PUBLIC bool_t bKXP84startRead()
{
    bool_t bOk = TRUE;

    uint8	read_address = 0x00;
    bOk &= bSMBusWrite(KPX84_ADDRESS, read_address, 0, NULL);

    return bOk;
}

PUBLIC uint8* i16KXP84readResult()
{
    bool_t bOk = TRUE;
    bOk &= bSMBusSequentialRead(KPX84_ADDRESS, 6, WorkBuffer);
    return WorkBuffer;
}

ソースコードの説明

KXP80のI2C仕様の12,13ページがメインです。
今回の手順はこの資料のSequence 4に相当します。
レジスタマップは資料後半にありますが、アドレス0から2バイト単位にX,Y,Zの加速度が格納されるようです。
http://akizukidenshi.com/download/KXP84-2050_Specifications_Rev2.pdf
まずMasterが通信開始の合図を出し、SAD+Wを出力します。
その後ACKを受けたMasterはRA(Register Address)を出す必要があるので、このアドレスを指定しておきます。
ここまでbSMBusWriteでやってくれます。
続いてi16KXP84readResultでSrとSAD+RをMasterが出してあげると、先の手順で指定したレジスタアドレスから順にデータを送ってきてくれます。
ACK/NACKの処理はbSMBusSequentialReadでやってくれます。

回路を作る

KXP84-2050からみたつなぎ方を言葉でメモしていきます。
VCCとGNDはそのままつなぎます。
MOTとFFは今は使わないので一旦放置。自由落下と過度加速度の検知に使うらしいです。
SCL/SCLKをTWE-LITEのSCLに
SDA/SDOをTWE-LITEのSDAに
ちなみにこの２つの線がI2Cなのでプルアップ必要かなと思いましたが、完成済みKXP84モジュール内部でプルアップされているので不要っぽいです。
モジュール内部のプルアップを抵抗の先はIO_VDDになっているので、こいつを持ちあげておきます。
CSは今回は使いませんが、プルアップしておきます。
RESETもプルアップ。あとでTWE-Liteからリセット制御できるようにしましょう。
ADDR0はプルダウンしておきます。この信号線でKXP84のI2C-SlaveAddressの最下位ビットが決まるとのこと。(つまりプルダウンした状態であればアドレスは0x18)

動作確認

ファームウェア書き込み後、PC側で端末(TeraTerm)を起動して設定をTWE-Liteにあわせておく。
起動させてターミナルで出力を確認

read finishX=[9270],Y=[AE38],Z=[8F30]
read finishX=[9270],Y=[AE38],Z=[8F30]

２回採取してみましたがよさげな値です。
センサを傾けると値が変わります。
アナログで測定した値とだいたい合っているので、一旦これで良しとします。

感想

今回とれた値をTWE-Liteのパケットで送信してあげて、受けた側(魔法陣モジュール側)で適当にLED制御すれば目的のものが作れそうです。
デバッグがソースコードと睨めっこになってしまうのがちょっと辛い。オシロかロジアナ欲しいです。
いっそ完成品買えば早いかもしれませんが、TWE-LiteからのI2Cの扱い方の勉強になったということで良しとします。
http://akizukidenshi.com/catalog/g/gM-09453/

最後にどうでもいいこと

TWE-Liteにファームウェアを書き込むとき、ライターに間違えたファイルを放り込むとこんなメッセージが出ます。
f:id:Catalina1344:20160201231155p:plain
ちょっぴりお茶目っ気があってクスっとしてしまいました。