研究

ブルースカイ・カタストロフィ

ブルースカイ・カタストロフィ(つよそう)

スクリーンショット 2018-05-05 17.10.41

↑は参考にあった Blue sky orbit in the Gavrilov-A. Shilnikov model でパラメータも同じもの.

以下は, \(\varepsilon\) をそのまま固定しながら, \(\mu\) を \([0.1, 1.1]\) の範囲で動かしたもの(0.02 間隔で 80000 回のルンゲクッタ.最初の 1000 回は捨てている).

OpenCVでカメラキャリブレーション

OpenCVを使ったカメラキャリブレーションをググると,まずCでの実装例が出てきて,次に比較的新しいpythonでの実装例が出てくるが,C++での実装例がいまいちヒットしない.

加えて,カメラキャリブレーションは基本的にあらかじめ用意したデータについて行うが,私は普段Webカメラを利用するので,データの収集とキャリブレーションを同時にできるといいな,というのがあった.

ので,作った.

リポジトリはこっち
https://github.com/Lait-au-Cafe/calibration

本当はできるだけシンプルにしたかったのだが,冒頭で述べたようにプリセットのデータがなくてもできるようにしつつ,でもプリセットのデータがあってもできる,という仕様にした結果少し複雑になってしまったように思う.

当該リポジトリを利用する際にはビルドにOpenCV3.xとpkg-configが必要なのに注意.

P行列の自由度とパラメータ

物の本によると,カメラによる物体の投影を表現する “P行列” は以下のようにあらわされる.

PMat1

なるほどな.

別の本では以下のように記述されている.

PMat2

まあそういうこともあるだろう.

ところで,P行列の自由度(DOF, Degree of Freedom)は一般に11とされる.
P行列は全部まとめると3×4行列になるが,定数倍しても意味が変わらないのでここから1を引いて11,という説明がされているのを見かける.

さて,上の二つの式をもう一度見てみよう.
回転行列Rと並進ベクトルtは二つの式で共通で,自由度はともに3であるから,一番左の行列だけで自由度は3+3=6.
真ん中の行列は焦点距離fがあって自由度1.
左の行列は,上の式では変数がa, s, c_x, c_yの4つで自由度4.
下の式では変数がδ_x, δ_y, c_x, c_yの4つで自由度4.

上の式も下の式も全部合わせて4+1+6=11自由度!
世界は今日も平和だなぁ…





ほんまか???

上の式と下の式で異なる点は以下の二点である.

  • 上の式ではx軸方向のスキューsが考慮されている.
  • 上の式ではアスペクト比aで表しているものを,下の式ではピクセルの物理的なサイズδ_x, δ_yで表している.

これを見て私が抱いた疑問が以下の2点.

  • 上の式でもaではなく物理的なサイズδ_x, δ_yで表せば自由度が12になる??
  • さらにy軸方向のスキューも考慮すれば自由度が13になる???

式にするとこんな感じ.

PMat3

しかし自由度12はギリギリ許容できても,3×4行列の自由度が13になることは逆立ちしてもあり得ない.
ではどれが正しくてどれがどう間違っているのか.

順に確かめていこう.

アスペクト比aとピクセルの物理的なサイズδ_x, δ_yについて



結論
“焦点距離f”と”アスペクト比a”のペアで表すか,”ピクセルの(相対的な)物理的サイズδ_x, δ_y”で表すのが冗長の無い表現である.

先程のP行列の式で,下の方の式では焦点距離fとピクセルの物理的なサイズδ_x, δ_yを用いてP行列を表現していた.
これは解釈上はいいのだが,自由度の観点ではこのうちの一つは冗長なので誤解を生じやすい.
変数はf, δ_x, δ_yの三つあるが,実際には自由度は2しかない.

式的な解釈としては,先ほどの式を少し変形して以下の形にするとわかる

PMat4

ここの式を見ると,1/δ_x’=f/δ_x, 1/δ_y’=f/δ_yと置き直しても普遍性を失わないことがわかる.
よって,スキューを考えない場合P行列の自由度は10.
あるいはf’=f/δ_x, a=δ_x/δ_yと置きなおすと,焦点距離とアスペクト比で表現できる.

f, δ_x, δ_yの変数のうち一つが不要である理由を言葉で表現するとどうなるか.
ここでもう一度焦点距離fの意味について問い直すと,これはカメラ座標系のスケールを正規化するためのパラメータとして解釈できる.
焦点距離fによる正規化を行った後,再度ピクセルの物理的なサイズδ_x, δ_yを用いてスケーリングを行うのだが,別に一度正規化を挟まずとも直接,カメラ座標系のスケールに対するピクセルの物理的なサイズδ_x’, δ_y’によってスケーリングしてしまえばそれでよい.
わざわざ一度fでスケーリングを行う必要はない.
(もちろん,処理に物理的な解釈を与えるという観点では重要な操作である).

もしスキューを考慮しないP行列で自由度が11になると説明する人間が居たら「ほんまか??」をぶつけよう.

y軸方向のスキューについて

結論
y軸方向のスキューは回転R,並進t,焦点距離f,アスペクト比aを取り直すことで消えるので考慮しなくてもよい.仮にy軸方向のスキューを明示的にP行列の式に組み込んでも自由度は変わらず11となる.

証明は以下.

PMat5

かなり雑に言葉で説明すると,カメラを90度回転させればx軸方向のスキューがy軸方向のスキューになる,という感じ.

画像作り終えてからsinの符号が逆なのに気づいた….
眠い…,適宜読み替えて読んでください….

まとめ

P行列の自由度は,
焦点距離f + アスペクト比a + スキューs + 画像中心の座標c_x, c_y
+ カメラの回転R(3自由度) + カメラの並進t(3自由度)= 11自由度.

以上.

OpenCL with OpenCV (OpenCV-CL)を使ってみた.-RGB編-

前回は取得したフレームをcvtColorでグレー画像に変換してからOpenCL側へ渡す,というまどろっこしいことをしていましたが,実際カラー画像を扱えないとお話にならないので何とか使えないかと試行錯誤.

CUDAの場合はuchar3は構造体として実装されているようで[要出典],uchar(=グレー画像)の場合とほとんど同様に扱うことができるが,OpenCLのuchar3は少し違った仕様になっているようで,ucharの場合をそのまま拡張しただけではうまくいかない.(CUDAのコードも適当に整理して上げたいですね).

ucharなどをスカラ型と呼ぶのに対し,uchar3などはベクタ型と呼ばれますが,このベクタ型の配列にアクセスする場合はvloadnやvstorenを使う.

参考サイト
How vector pointers work in openCL
Vector Data Load and Store Functions

それを考慮して実際にRGB画像を扱って各要素にアクセスした例が以下.

やっていることは至極単純なグレースケール化のみ.
今回もよくわかっていない点として,widthがいつでもUMatのpitchと等しくなるのか,という点.
例えば,CUDAの場合はx+y * widthとやるとうまくいかなくて,x+y * pitchのようにしなければいけない.つまり要素数はwidth分だが,各行当たりのメモリはもう少し余裕をもって確保されている.このpitchの値はcudaMallocPitchでメモリをアロケートする際に一緒にもらえる.
今回のOpenCV-CLの場合はOpenCV側の提供するUMatを使っているため問題ないのだろうか.
一応UMatもプロパティとしてstepやoffsetを持っているため,これを考慮せずにindexingしても大丈夫なんだろうかいやでも今のところうまくいってるしなあ,という感じ.
つまるところ,まだまだ検証の足りないコードなのでindexing周りで不具合が生じる可能性がある.

ちなみにindexingまわりで不具合があると出力画像はこんな風になりがち.
DXsiYnIVwAAS8bV

今日は以上.

OpenCL with OpenCV (OpenCV-CL)を使ってみた.

ここに “””つらみ””” があります.

image

環境を統一しようにもなかなか険しいものがある.(NvidiaのOpenCLはCUDAのドライバ入れれば動く…?)
世界に試されている.
でも今出先なので(=Geforceを積んでいるデスクトップ環境が使えないので),ノートPC上でできる環境を適当に決めて実装しないと無限に「グレブナ基底と代数多様体入門Ⅰ」を読んでしまう.無限に某ーロッパ企画のゲームをプレイする某員長の動画を見てしまう.

色々と考えた結果,OpenCVの環境作ればOpenCLが使えるらしいことを知って,じゃあこのOpenCL使ったらいいか,という気持ちになった.
OpenCVは全体的にドキュメンテーションに難があってほとんど参考にならないが,神先駆者様がいたのでこれをほとんどコピペでとりあえずサンプルコードを動かしてみた.

negaposi

とりあえずサンプルは動いたが今後手を加えていくとうまくいかなくなるかもしれないのでしばらく様子見.
カーネルコードをロードしてコンパイルするあたりでRustのシャドーイングを使いたくなったが残念ながらこれはC++なのでそうは問屋が卸さない.
世界中のスクリプトがすべてRustに置き換わらないかな.
せめて世界中の言語がすべて強い静的型付き言語にならないかな.

RustからC++を呼び出すこともなんとなく考えたが,参考サイトに
“C++ exceptions can mostly be caught in D, but that’s not a thing in Rust at all. “
とあって,確かになぁ…という気持ちになった.
ただでさえセーフな言語からアンセーフな言語を呼び出すのは厳しいのに,エラーハンドリングのやり方などの言語仕様も異なるとなると,考えることが増えそう(曖昧).

コードの印象としては,カーネルコードのロード,コンパイル,呼び出しの流れはPyOpenCLのほうが若干すっきりしていた気がする.(以前書いたPyOpenCLではホストコード中にカーネルコードを書いていたのに対し,今回は別ファイルにカーネルコードを書いているため,ファイルを読み込む処理が追加されているため煩雑に感じるのかもしれない).
カーネルコードの引数指定の部分は知らないと絶対にわからない仕様(詳しくは参考サイト様参照)なのにドキュメンテーションには記載が見当たらなかったので,このあたり敷居が高いなと感じた.

今回のコードで少し不思議に感じたのは,フレームのグレースケール画像(gray)が更新されたときにデバイスのバッファ(d_frame)へデータを送らなくてもちゃんと更新が反映された点.UMatは今回の例のように直接imshowに渡せるようだし,情報の受け渡しはOpenCV側がやってくれているのだろう.
つまり,cv::UMat d_frame = gray.getUMat(…);が実行された時点で,grayへの書き込みが自動的にd_frameに反映されるようにお膳立てしてくれているのかもしれない.(というかそもそも実体がGPU上に移るのかもしれない).
確信はないけど.

とりあえず今日はここまで.

参考サイト様
OpenCV 3.0.0での独自カーネルOpenCL
UMatの内部処理(OpenCL編)
Any updates on calling C++ from Rust

unityで外部テキストをC#コードとして読み込む

こんにちは。mo-takusanです。

今まではほとんど記事を書いてこなかったのですが、さすがにそろそろ投稿していこうと思います。

unityで外部コードを読み込むには…

最近とある事情でコンパイル後に外部ソースコードを読み込みたい、という場面がありました。

「あああ、Python書くかぁ」

となっていたのですが、念のため調べると、C#はコンパイル言語ですが、ランタイムでコンパイルができるようです。C#神!!やっぱり.NET最強!!!

ただし、untiyとの相性が悪いものが多く、四苦八苦してしたので、紆余曲折の様子を記事にしていきます。

なお、本記事で使用したunityのバージョンは2017.1.1f1です。

まずはできた方法を、、

まあ、とにかく知りたいのは実際にうまくいったやり方でしょう(需要があれば)から、その方法から話していきます。

なんだかんだとググっていくうちに以下のGitHubにたどり着きました。

https://github.com/aeroson/mcs-ICodeCompiler

このサイトを見たとき、ついに来た!最終コミットもそれなりに新しいしこれはいけるやろ!と思ったのですが、ここからが長かった…(unityのバージョンを合わせればよかったとかそういうツッコミはなしでお願いします、、、)

2018-02-07_23h38_41

早速クローンしてサンプルプロジェクトを開いてみると…

2018-02-07_23h42_07

地獄絵図…。unityはICodeCompilerやCompilerResultsが存在していないと主張してくるのですが、

https://msdn.microsoft.com/ja-jp/library/system.codedom.compiler.icodecompiler(v=vs.110).aspx

https://msdn.microsoft.com/ja-jp/library/system.codedom.compiler.compilerresults(v=vs.110).aspx

のように、System.CodeDom.Compiler名前空間には確かに存在しています。実際コードの方でもusingされており問題はなさそうです。やけくそで丁寧に書き直してみても全然ダメ…

2018-02-07_23h45_47

うーーーん。

これについては解決方法が全く分かりませんでした。ので、ダメもとで一緒にクローンしてきた「non Unity3D usage」のコードを新規プロジェクトに入れてみました。dllはPluginsに入れることを忘れずに…

2018-02-08_00h03_29

おや?エラーが出ない!
やった…やったよ、、おれ、ついにやってのけたよ…
そうするとデモが動かなかった理由がいよいよ謎ですね。

それはともかくとして、あとは適切にファイルを読み込むコードを書いて、、、

 private static Assembly Compile(IEnumerable&lt;FileInfo&gt; files)<br />
 {<br />
     var domain = AppDomain.CurrentDomain;<br />
     var references = domain.GetAssemblies().Select(a =&gt; a.Location).ToArray();<br />
     var options = new CompilerParameters<br />
     {<br />
         GenerateExecutable = false,<br />
         GenerateInMemory = true,<br />
     };<br />
     options.ReferencedAssemblies.AddRange(references);<br />
     var compiler = new CodeCompiler();<br />
     var result = compiler.CompileAssemblyFromFileBatch(options, files.Select(f =&gt; f.FullName).ToArray());<br />
     return result.Errors.Count &gt; 0 ? null : result.CompiledAssembly;<br />
 }

完成しました!!!
なお、こちらのプロジェクトは諸事情で公開していないのであしからず。

ここからはできなかったことを書くぞ!

と言いましたが、疲れたので失敗の様子はダイジェストでどうぞ

System.CodeDom.Compilerを素で使う

なぜかコンパイルすると動かない(unityではよくある泣)

AssetStoreからダウンロード

https://assetstore.unity.com/packages/tools/cs-script-for-unity-23510

このページにいい感じのアセットがあったのでダウンロード。が、だめ。(古かった)

Roslynを使う

unityが.NET4.6に対応したことを利用して、次のサイトを参考にして組んでみよう!

https://www.gaprot.jp/pickup/tips/roslyn

unityが反応してくれない!おい!!nugetからdll落としてPluginsに入れてもダメでした。(むしろそれ以外を知りませんが)

ローレンツアトラクタとカオス

こんにちは.TRSasasusu です.

記事を更新していないだけで,KCS は活動しています.新歓の準備も進んでいるようです.

さて,自分は生命情報学科に所属していますが,この前のレポートにてカオスで遊んだので,今回はその話をします.

スクリーンショット 2018-02-01 15.46.11

これはローレンツ方程式(↓の式)のアトラクタです.この軌道に沿って点が移動していきます.


\begin{align}
\frac{dx}{dt} &= -\sigma x + \sigma y \\
\frac{dy}{dt} &= -xz + rx -y \\
\frac{dz}{dt} &= xy – bz
\end{align}

特にこのようなフラクタル構造を持つアトラクタはストレンジアトラクタと呼ばれ,カオスを生み出します.蝶の羽みたいで美しいですね.

そして,パラメータ(今回は \(\sigma,\ r,\ b\))を変化させると,カオスになったりならなかったりします.これを表すのが分岐図です.

スクリーンショット 2018-02-01 16.07.43

なんかすごいですね.そして,カオスを定量的に表すのがリアプノフ指数(多次元でまとめたものがリアプノフスペクトラム)です.力学系を表すヤコビ行列を時系列順に並べて QR 分解して(時刻 \(0\) では時刻 \(0\) のヤコビ行列のみを,それ以降ではヤコビ行列に前時刻の \(Q\) を右から掛けたものをそれぞれ QR 分解する),各時刻の \(R\) のそれぞれの対角要素の \(\log\) をとって全部足し合わせて \(N\) で割ると,以下のようになります.

スクリーンショット 2018-02-02 1.06.55

1 つでも正のリアプノフ指数があればカオスとなるようなので,この結果は分岐図のものと一致していますね.

最後に,ストレンジアトラクタの動画を載せて終わりにします(動画はレポートには載せられない).

インスタ度判定Twitterアカウント@insta_san

今年はAI班で三田祭に向けてインスタっぽい画像を判定するアプリを作成しました。その名もインスタさんです。

インスタさんに画像を渡すと・・・

Screenshot from 2017-11-23 00-52-32

このように画像のインスタっぽさと、インスタっぽい部分が表示されます!
(三田祭期間中のみの運用予定なので、記事を閲覧しているときには動かないかもしれません、ごめんなさい!)

AIは@sesenosannko、Twitterの入出力は@mt_caretが作成しました。

しくみ

インスタさんのアイデアはとても簡単です。Instagramからたくさんの画像といいね数を集めてきて、いいね数が多い画像をインスタっぽい画像として学習させます。

ここで考えなくてはいけないことは、「いいね数が多いという基準」と「学習モデル」です

  • データセット

画像からインスタっぽさを学習するのは難しい課題なので、今回は単純に「インスタっぽい」画像と「インスタっぽくない」画像の二つに分類することにしました。単純にいいね数が多い画像と言ってもユーザーによってフォロワー数が違いますし、投稿した時点でのフォロワー数がわからないので、いいね数をそのまま使うことは難しいでしょう。

調査の結果、多くの人のアカウントでは投稿するたびにいいね数が徐々に増えていく傾向にあることが分かりました。しかし、人によって増加のタイミングや速度は様々で何かの関数で近似するのは難しそうだったので、単純に前後20投稿の平均いいね数よりもいいね数が多いか少ないかによって「インスタっぽい」画像と「インスタっぽくない」画像を分けることにしました

  • 学習モデル

近年は画像の分類にはCNNというディープラーニングのモデルが使われることが多いです。インスタさんもCNNを用いて学習しています。

今回は画像数が少ないこともあり、大量の画像で分類を学習した結果を活用して学習を行なっています。このような手法は転移学習と言われています。使用したモデルの詳しい説明は以下のサイトを参照してください。

https://www.tensorflow.org/tutorials/image_retraining

画像がインスタっぽいかという問題は人間でもはっきりと分けられないので、あまり良い学習ができることは期待できません。特に今回は学習セットの分け方も非常に単純で信頼性の低いものとなってしまっています。しかし、転移学習によって画像から意味のある情報を取り出して学習を行なっているため、何かを学習しているということは見て取れます。

LIME

インスタっぽい部分を出力するときに使用しているのはLocal Interpretable Model-Agnostic Explanations(LIME)という手法です。

簡単にいうと、入力画像の一部を隠した画像をたくさん作り、インスタ度が高かったときに隠されていなかった部分のインスタ度が高いと考える手法です。詳しくは以下のサイトを参照してください。

https://homes.cs.washington.edu/~marcotcr/blog/lime/

StackTraceを実装する話

* この記事では、MSVCのinlineアセンブリ構文を使う。
** C++はあまり詳しくないのと、深夜テンションで思いついた実装なので、もっと速い・簡単なやり方があればぜひ教えてほしい。

Unity3Dをやっている人は分かるが、簡単に実例を説明しよう。例えば、描画関数の中でしか使えない関数がある。

void OnDraw () {
    Game::DrawCube(...);
}

ここでGame::DrawCubeはOnDrawおよびOnDrawで呼ばれた関数からの実行出なければいけない。では、エラーチェックとしてOnDrawから呼ばれたかを確認したい。

じゃあ、これをどう実装すればいいのかという話。もちろんStackUnwindとかのライブラリ使ってもいいけど、そもそもデバッグではなくて実際のリリースでも使いたいのと、関数の名前より関数のIDを見るのと、個人的に人の実装を使いたくない癖(ここ重要)がある。

  • StackFrameについて

さて、C++(別にC++じゃなくても)のDisassemblyに着目する。関数が呼ばれたときにstackをみると

address   value
------------------------
EBP-4     return address
EBP       old EBP
EBP+4     funcvar1
EBP+8     funcvar2
...

になっている。これを関数のStackFrameと言う。

では、EBPはこの関数のヘッドを指しているが、その値がなんと呼ばれた関数のヘッドのアドレスが入っている。つまり、

mov EBP [EBP]

をやれば、StackTraceをすることができる(EBPは破壊されるのでこのまま使わないが)。

  • アドレスに着目

関数というのは、あくまで下のようにStackにある変数を操ったりする命令の集合である。例えば、以下ではDrawCubeを呼ぶOnDraw関数の例である。

        OnDraw:
0x123     instruction 1
0x125     instruction 2
...
0x200     call Game::DrawCube (0x502)
0x204     ret

(0x502はGame::DrawCubeのラベルが存在する位置である。例えばOnDrawだと0x123となる)

ここで、Game::DrawCubeが呼ばれた時のStackFrameを見ると、return addressが0x204となる。つまり、この関数が終わったらEIP(instruction pointer)がどこを指せばいいかを教えている。

では、この情報をどう使えばいいかとなるが、まず「呼ばれた関数」の定義を改める。「この関数はあの関数であるか」というのは、厳密に言うと「この関数のreturn addressは一緒なのか」となる。なぜかというと関数の名前より、システムが描画したいときにしか呼ばない「一ヵ所」からの実行だからである。

  • 実装

まず、OnDraw関数を呼ぶところは一定なので、下のようにあらかじめ登録することができる。なお、ラベルはアプリで共通されるため、__COUNTER__などでユニークなラベルにする必要がある。

int funcLoc;
__asm{
    uniquelabel1:
    mov EAX, offset uniquelabel1
    mov funcLoc, EAX
}

なお、このコードは関数を呼んだ直後に入れる。

次に、DrawCubeの中で、親(の親の親の…)のreturn addressがfuncLocであるかを確認すればいい。

偽コード
for EBP=*EBP
while *(EBP+4) != funcLoc ;

実際のコード

__asm{
    mov EBX, EBP
  loop:
    mov EBX [EBX]
    mov EAX, [EBX-4]
    cmp EAX, funcLoc
    jne loop
}
  • Access Violation防止

上のループは、funcLocが見つからないと永遠にループしてAccessViolationエラーが起こる。これを防止するにはどこかで止める必要がある。まず、mainスレッドからの関数と考えれば、
1. assert (std::this_thread::get_id() == mainThreadId);
2. asmのなかに

cmp EAX, mainThreadFuncLoc
je end

ただし、mainThreadFuncLocはmain()を呼んだ__CRT_Startupなんちゃらの場所である。

  • まとめ

アセンブリをいじることで、関数の親をトレースすることができた。この方法を使って、この関数を呼んだ親を区別することも可能となる。なお、デバッガでRTCが入ると、関数が実行された後にEBPの確認instructionが入るので、アドレス登録が動かなくなる場合がある。解決方法考え中(とりあえずオフにした)。もちろんリリースでは影響されない。

Winapiでプロセス間のリアルタイム映像転送

どうも、チョコです。久しぶりの進捗がありました。

今回は、OpenGLの2つのアプリ間で通信をしないといけないことになりました。ただし、リアルタイムの生画像転送なので、100mbpsくらいの速度が最低ないと困る。

ポートとかパイプあけて通信すればよいではという考えもあったが、アプリAがアプリBの好きなデータを取れるようにしたい(そういうプロジェクトだから)。なので、いちいちデータの名前送ってデータをもらうとすごく面倒だし、なにより遅い。

そこで、WinAPIのReadProcessMemoryです。要は、アプリAがアプリBに親権限を持っていれば、なんでもしていいよとのこと。よさ。
WriteProcessMemory(HANDLE processHandle, void* source, void* destination, ulong numbytes, ulong* numbytesread)

*sourceはアプリBでのポインタ値で、destinationはそのポインタの指しているデータをアプリAにコピーする場所。destinationを先に初期化しないといけない(それはそう/1時間消耗した。)。
ん、読みたいところのポインタがないとできないね(それはそう)。ので、まずはポインタ集というStructを作ります。
struct PointerList {
uint hasDataLoc, pixelsLoc, pixelCountLoc, displayWLoc, displayHLoc, okLoc;
} pointers;

*win32なので、ポインタは32ビットunsigned intに変換できます。
なので、最初にこのStructを読み取れば、すべてのポインタがわかるね。わーい。

実際にやったものを見せよう。

わーい。(ただこれを見せたかった)

では。