2013年12月10日火曜日

Ceres Solverの使用例(放物線の最小二乗法)

Ceres SolverはGoogleが開発した非線形最適化のためのC++ライブラリ。 使い方は本家サイトにもチュートリアルがあるし、サンプルコードもある。 でも、どれも単純すぎたり、難しすぎたり中途半端なので、 データが復数ある最小二乗法のサンプルコードを載せる。

\[ y=a x^2 \]

この関数をフィッティングする場合、データとして$(x_i, y_i)\ (i=0,1,\cdots,n-1)$が与えれるので、 求めるパラメータは、$a$のみ。ここのデータの誤差を \[ e_i = y_i - a*x_i^2 \] とすると、単純な最小二乗法では、 \[ E=\sum_i e_i^2 \] が最小化される。 ここでは、データ$(x_i, y_i)$として次のファイル「data.txt」が与えられるとする。

0 0
1 1.2
2 3.9
-1 0.9
-2 4.1

サンプルコード

まず、ヘッダー部分は単純にceresのヘッダをインクルードすれば良い。 namespaceはライブラリのexamplesを真似して追加した。 OpenCVのヘッダはデータ入力に使うだけなので本来必要ない。

#include "ceres/ceres.h"
#include "opencv2/opencv.hpp"
 
using ceres::AutoDiffCostFunction;
using ceres::CostFunction;
using ceres::Problem;
using ceres::Solver;
using ceres::Solve;

次に、誤差関数の定義部分。 ErrorFunctionクラスのメンバ関数として定義する必要がある。 operator()関数をテンプレート関数として ここに書くのは$e_i$の部分だけ。 パラメータ$a$も残差$e_i$も配列で定義するみたい。 汎用性を考えてType型という抽象型になっているだけで今回はdoubleと思えば良い。

class ErrorFunction{
private:
    //1組のデータ
    const double x;
    const double y;
 
public:
    //データを1つずつ入力する際このコンストラクタが呼ばれる。
    ErrorFunction(double x, double y) : x(x), y(y) {}
 
    //e_iの定義。
    template <typename Type> bool operator()(const Type* const a, Type* residual) const {

        //e_i = y_i - a*x_i*x_i
        residual[0] = Type(y)  - a[0]*Type(x)*Type(x);
        return true;
    }
};

誤差関数が定義できれば、あとはデータを与えて最適化するだけ。 データの入力は面倒なので 「cv::Mat_をテキストファイルに保存・読み込み 」 で紹介した方法を使う。

void main(int argc, char *argv[]){

    //最適化するパラメータ。初期値を入力しておく。
    double a = 0.0;
 
    //データの入力
    cv::Mat_<double> data;
    readTxt("data.txt", data);
 
    //データの入力
    Problem problem;
    for (int i = 0; i < data.rows; ++i) {

        //x_i,y_iを1つずつ入力する
        problem.AddResidualBlock(
            new AutoDiffCostFunction<ErrorFunction, 1, 1>
            (
                //ここでデータを入れている
                new ErrorFunction(data(i, 0), data(i, 1))
            ),
            NULL,    //Loss function
            &a
        );
        //std::cout << data(i, 0) << ' ' << data(i, 1) << std::endl;
    }
 
    //ソルバの準備
    Solver::Options options;

    //CXSparseをインストールしていないときはDENSE_*にしないとエラーになる。
    options.linear_solver_type = ceres::DENSE_NORMAL_CHOLESKY;   
 
    //最適化
    Solver::Summary summary;
    Solve(options, &problem, &summary);
 
    std::cout << "Fit complete:   a: " << a << "\n";
    return;
}

参考

Comparison of head-free gaze trackers

This page summerizes head-free gaze tracking hardwares.

SteelSeries Sentry Gaming Eye Tracker

https://steelseries.com/gaming-controllers/sentry-gaming-eye-tracker

  • Price:~$250
  • Acquisition rate:30 or 60 Hz
  • Angular resolution:0.5 - 1 deg (RMS)
  • Latency: < 20ms (at 60Hz)
  • FOV: 40cm x 30cm(at 30Hz, depth 65cm)
  • Depth range: 45cm – 75cm
  • Size: 20 x 1.9 x 1.9 cm
  • Weight: 70g

EyeTribe

https://theeyetribe.com/

  • Price:~$100
  • Acquisition rate:30 or 60 Hz
  • Angular resolution:0.5 - 1 deg (RMS)
  • Latency: < 20ms (at 60Hz)
  • FOV: 40cm x 30cm(at 30Hz, depth 65cm)
  • Depth range: 45cm – 75cm
  • Size: 20 x 1.9 x 1.9 cm
  • Weight: 70g

Gazepoint GP3 Desktop Eye Tracker

http://gazept.com/

  • Price: ~$500
  • Rate: 60Hz
  • Angular resolution:0.5-1 deg (RMS)
  • Latency: 50 – 70ms
  • FOV:25cm x 11cm (at 65cm)
  • Depth range: ±15 cm
  • Size: 20 x 1.9 x 1.9 cm (W x H x D)
  • Weight: 70g

S2 Eye Tracker

http://mirametrix.com/products/eye-tracker/

  • Price: ~$500k
  • Rate: 60 Hz
  • Angular resolution: 0.5 deg (RMS)
  • Latency: ? (0.3 deg)
  • FOV: 25 x 11 cm
  • Depth range:30 cm
  • Size: 35 x 4 x 3 cm (W x H x D)
  • Weight: 300 kg

Tobii X2-30 Eye Tracker Compact Edition

http://www.tobii.com/de/eye-tracking-research/global/products/hardware/tobii-x2-30-eye-tracker/

  • Price: ~$1M
  • Rate: 30 Hz
  • Angular resolution: 0.4 - 0.5 deg (RMS)
  • Latency: 50 - 70ms
  • FOV: 50 x 36 (at 70 cm)
  • Depth range: 40 - 90 cm
  • Size: 18.4 x 2.8 x 2.3 cm
  • Weight: 200g

Tobii X2-60 Eye Tracker Compact Edition

http://www.tobii.com/zh-CN/eye-tracking-research/global/products/hardware/tobii-x2-60-eye-tracker/

  • Price: ~$1M
  • Rate: 60 Hz
  • Angular resolution: 0.4 - 0.5 deg (RMS)
  • Latency: 35 ms
  • FOV: 50 x 36 (at 70 cm)
  • Depth range: 40 - 90 cm
  • Size: 18.4 x 2.8 x 2.3 cm
  • Weight: 200g

2013年12月7日土曜日

Ceres Solverのコンパイル方法 on Windows

Ceres Solverとは?

Ceres SolverはGoogleが開発した非線形最適化のためのC++ライブラリ。 何ができるのかというと、定義した誤差関数を最小化するパラメータを求めることができる。 非線形最適化のための定番アルゴリズムLevenberg–Marquardt法(LM法)などは実装済みで利用者は殆ど意識すること無く使用することができる。 何よりも、Automatic Differentiationという機能がすごい。 誤差関数の微分を数値的にではなくて、解析的に(数式的に)自動的に解いてLM法に組み込んでくれる。 つまり誤差関数そのものを実装するだけで良い。

誤差関数の微分の計算がなくなっただけでそんなに楽になるのだろうか? 例えば、ある誤差関数$f(x,y)$を作成し、これを最小化したいとする。 効率よく計算するためには、この関数の解析的な偏導関数$\frac{\partial f(x,y)}{\partial x}$と$\frac{\partial f(x,y)}{\partial y}$が必要になる。 この場合変数は$x$と$y$の2つだけど、変数が増えると必要な偏導関数の数が増える。 画像自体を推定する場合は、全ての画素の輝度がパラメータになる。 for文などで自動化出来ても$\sin$、$\cos$関数などが含まれている場合は、個々の偏導関数自体が非常に複雜になる場合もある。 誤差関数が複雜な場合、手で計算するのには限界がある。 Mathmaticaのような数式処理ソフトを使う方法もあるが、実装もなかなか面倒だった。

こうした実装をC言語でかつ安直な方法でやろうとすると、 読みづらいマクロと関数ポインタ、グローバル変数を使いまくる実装になってしまう。 Ceres SolverはC++の継承と関数オブジェクトを使って非常に美しく、使いやすい形でまとまっている。 現在のところ日本語の解説は、LinuxのみだったのでWindowsでの動作実績をまとめる。

インストール手順

  1. 必要なファイルライブラリをダウンロード、ビルド
    1. eigen-eigen-ffa86ffb5570.zipをダウンロード(3.1以上)する。ビルドの必要はない。
    2. glog-0.3.3.tar.gzをダウンロードする。
    3. 解凍し、oogle-glog.slnを叩いてビルドする。バッチビルドで全てを選んでビルドする。エラーが出るがビルドは成功。
    4. gflags-2.0-no-svn-files.tar.gzをダウンロード。問題無くビルド成功。
    5. 解凍してgflags.slnの読み取り専用設定を解除。
    6. gflags.slnを叩いてビルドする。
  2. ceresのダウンロード、フォルダ作成
    1. ceres-solver-1.8.0.tar.gzをダウンロードする。
    2. ceres\ceres-install、ceres\ceres-buildフォルダを作成する。
  3. cmakeによるVisual Studioのソリューションファイル作成
    1. cmakeを起動する。
    2. 「Where is the source code」にceres/ceres-solver-1.8.0設定する。
    3. 「Where to build the binaries」にC:/library/ceres/ceres-buildを設定する。
    4. 一度configureする。エラーが出て止まる。
    5. GroupedとAdvancedにチェックを入れると全てのオプションが見えるようになる。
    6. CMAKE_INSTALL_PREFIXにceres\ceres-installへの絶対パスを入れる。
    7. EIGEN_INCLUDE_DIRにeigen3.2への絶対パスを入れる。
    8. 再度configureする。エラーが出て止まる。
    9. GLOG_INCLUDE_DIRにglog-0.3.3\src\windowsへの絶対パスを入れる。
    10. GLOG_LIBRARYにglog-0.3.\Release3\libglog.libへの絶対パスを入れる。
    11. GFLAGS_INCLUDE_DIRにgflags-2.0\src\windowsへの絶対パスを入れる。
    12. GFLAGS_LIBRARYにgflags-2.0\Release\libgflags.libへの絶対パスを入れる。
    13. CX_SPARSE_INCLUDE_DIRにCXSparse\Includeの絶対パスを入れる。
    14. CX_SPARSE_LIBにCXSparse\Lib\CXSparse.libの絶対パスを入れる。
    15. 再度configureする。
    16. generateするとceres\ceres-buildにCERES.slnが生成されている。
  4. CERES.slnを叩いてビルドする。バッチビルドでALL_BUILDとINSTALLのチェックを全て付けてビルドする。

テスト環境

  • Visual Studio 2010 sp1
  • Windows 7 Pro

参考

2013年11月28日木曜日

Gnuplotによるカメラ画像のPSF推定

カメラ画像のボケぐあいであるPSF(Point Spread Function、点拡がり関数)を求める1つの手法をまとめる。 Gnuplotによる区分関数の最小二乗法 を用いるとプログラミングの必要なく非常に簡単にPSFを求めることができる。 ここで仮定するPSFは、次の式で表現されるピルボックス関数(Pillbox function)である。

PSFとしてピルボックス関数を仮定する

\[ R(x, y, r)= \left\{ \begin{array}{cl} \frac{1}{\pi r^2} & x^2+y^2 \leq r^2\\ 0 & {\rm otherwise} \end{array} \right. \]

階段関数のテストパターン

テストパターンとして次の図のような白黒2色のパターンを撮影する。 これは水平方向に見れば階段関数(step function, staircase function)と見なすことができる。

\[ S(x,y)= \left\{ \begin{array}{cl} 1 & (x \leq 0)\\ 0 & (x < 0) \end{array} \right. \]

観測パターンとそのモデル

上記のパターンがピンぼけて下図(原画の4倍スケール)のようなパターンが観測される。

ピルボックス関数とテストパターンの畳み込み積分を考えると、 ピルボックス関数の一部がテストパターンの黒い領域に重なり、その残りの部分が白い領域に重なる瞬間がある。 この2つ重なり領域のピルボックス関数の体積比が画像の輝度として観測されていると考えることができる。 こう考えると観測画像の水平方向の輝度は下式で表されることが知られている。 $-1 \leq t \leq 1$の範囲がピルボックス関数とテストパターンの黒白境界が重なる範囲である。

\[ f(x)= \left\{ \begin{array}{ll} I_{\rm min} & t<-1 \\ g(t)(I_{\rm max} - I_{\rm min}) + I_{\rm min} & -1 \leq t \leq 1 \\ I_{\rm max} & t>1 \end{array} \right. \] \[ g(t) = \frac{1}{2} + \frac{1}{\pi}(t \sqrt{1-t^2}+\sin^{-1}t) \] \[ t = \frac{x - c}{r} \]

PSFのパラメータ推定

上記関数のパラメータは何か? パラメータは白い領域の輝度$I_{\rm max}$、黒い領域の輝度$I_{\rm min}$、中心位置$c$、ボケ幅$r$だけど、 PSFとして意味のあるパラメータは$r$だけである。 4つのパラメータを求めてそのうちボケ幅$r$だけがピルボックス関数の形を決めるパラメータとなる。

一方データとしては、画像の水平方向に走査したときの輝度値を用いる。 垂直方向に輝度値を平均すれば画像のノイズの影響を軽減することができる。 このデータに対して先ほどの関数をGnuplotを用いてフィッティングする。 フィッティングされた関数とデータの例を下に示す。

reset

# 関数の定義
t(x,c,r)=(x-c)/r
g(x,c,r)=0.5+(t(x,c,r)*sqrt(1-t(x,c,r)**2)+asin(t(x,c,r)))/pi
f(x,c,r,imax,imin)=(t(x,c,r)<-1) ? imin :((t(x,c,r)>1) ? imax : g(x,c,r)*(imax-imin)+imin)
 
# 初期値の設定
c=50
r=25
imax=95
imin=5
 
# フィッティング
fit f(x,c,r,imax,imin) 'sample.txt' via c,r,imax,imin
 
# データと関数描画
set xr[-5:105]
plot 'sample.txt' with points pt 3
replot f(x,c,r,imax,imin) linewidth 3

Gnuplotの出力は以下の通り。 以下の結果からボケ幅は$r=21.298$画素であることが分かる。

 Iteration 0
 WSSR        : 587.799           delta(WSSR)/WSSR   : 0
 delta(WSSR) : 0                 limit for stopping : 1e-005
 lambda   : 0.859572
 
initial set of free parameter values
 
c               = 50
r               = 25
imax            = 95
imin            = 5
/
 
 Iteration 1
 WSSR        : 68.9751           delta(WSSR)/WSSR   : -7.5219
 delta(WSSR) : -518.824          limit for stopping : 1e-005
 lambda   : 0.0859572
 
resultant parameter values
 
c               = 49.3057
r               = 21.7139
imax            = 94.3184
imin            = 6.75013
/

   :
  (省略)
   :
 
 Iteration 5
 WSSR        : 63.5581           delta(WSSR)/WSSR   : -8.02681e-007
 delta(WSSR) : -5.10168e-005     limit for stopping : 1e-005
 lambda   : 8.59572e-006
 
resultant parameter values
 
c               = 49.278
r               = 21.298
imax            = 94.2626
imin            = 6.80794
 
After 5 iterations the fit converged.
final sum of squares of residuals : 63.5581
rel. change during last iteration : -8.02681e-007
 
degrees of freedom    (FIT_NDF)                        : 96
rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 0.813673
variance of residuals (reduced chisquare) = WSSR/ndf   : 0.662063
 
Final set of parameters            Asymptotic Standard Error
=======================            ==========================
 
c               = 49.278           +/- 0.07501      (0.1522%)
r               = 21.298           +/- 0.1477       (0.6935%)
imax            = 94.2626          +/- 0.1459       (0.1548%)
imin            = 6.80794          +/- 0.1517       (2.228%)
 
correlation matrix of the fit parameters:
 
               c      r      imax   imin   
c               1.000 
r              -0.013  1.000 
imax            0.432  0.325  1.000 
imin            0.449 -0.338 -0.013  1.000 

参考

  • Gnuplotによる区分関数の最小二乗法
    http://buaiso.blogspot.jp/2013/11/gnuplot.html
  • B. Okumura, M. Kanbara, and N. Yokoya: "Augmented reality based on estimation of defocusing and motion blurring from captured images", Proc. IEEE and ACM Int. Sympo. on Mixed Augmented Reality (ISMAR 06), pp. 219-225, 2006.

2013年11月15日金曜日

左翼先生の栄えた時代

自分だけの稀な体験なのか、それとも全国的に似たような傾向があるのか全く分からないけど、 自分が義務教育時代に経験した、というより違和感を感じた、出来事をここで紹介したい。

はっきり言って高校に入るまで心から尊敬できる先生は居なかった。 その最大の原因は、先生があまり勉強していないことが、生徒である自分にもなんとなく分かったから。 専門の教科ですら、質問しに行けばろくに答えられずにごまかす先生も多かった。 生徒は子供だから大丈夫だろうというのは全くの間違いで、その時はそういうものだと思っても必ず覚えていて、 大人になってやっぱりあの先生があの簡単な質問に答えられなかったのは、おかしいじゃないかと気づくもの。 自分の場合も、そういえば、今考えてみたら、当時の先生は変な人が多かったなという印象をもつ。 今思えば、今の自分の知識からすれば、当時自分が感じていた以上にその先生方は変な先生だった。 どう変かと言えば極端に左に偏った思想をもち、ろくに勉強もせずに自分の思い込みを、 めちゃくちゃな論理で押し付けてくる。そんな先生があちこちに居た時代だった(今は知らない)。

めちゃくちゃな論理

小学校時代の遠足や運動会の時だけお弁当を持って行くことになっていた。 ただ、なぜかお弁当の中身に関して変なルールがあった。

お弁当はおにぎり2個と漬物以外持ってきてはならない。

この理由は今の人には全く分からないだろう。 要するに貧乏な人もいるからその格差が目立ってはいけないので、 全員貧乏な飯にすべしということだった。 いかにも平等が大好きな左翼先生が考えそうなことだ。 当時の自分は、お弁当も持ってこれないような貧乏がいるかよと思った。 確かに芸人のダウンタウンの子供時代の話を聞けば、塩むすびとみょうがしか持ってこない生徒も居たとか。 それに、共働きでお弁当を作ってくれない親もいるかもしれない。 そういう配慮かもしれないけど、全員を最下層に合わせるという発想が貧困というか、だったらパンでもいいじゃないかと思ってしまう。 子供の頃は人との違いを非常に気にして、ちょっと変わった人がいたらいじめてしまうものである。 でも、それは何も貧乏かそうでないかだけじゃなくて、体型だとか、顔だとか、名前ですら。 ありとあらゆるものに関してである。 決して、体型が分からないようにダボダボの服を着ましょうとか、名前はやめて優劣の無い記号にしましょう、とは言わない。 にも関わらず、弁当だけが平等の象徴かのように取り上げてくだらないルールを作るのは、 平等が大好きなあまり何が大事なのかということを考えなくなっている証拠じゃないかな。 この弁当ルールに関して、「育ち盛りの子供に米+漬物のみ食事をさせるなんてありえない」と言い、 大きなおにぎりの中に大量の具材を入れてくれた、母親の方がよっぽどまっとうな考え方のように思える。

中学校時代にも無意味な校則が存在した。 今でも納得行かないものは、冬場の手袋とウィンドブレーカーに関する校則である。 なぜか、手袋を着用できるようになる日が決まっていた。 しかも、かなり寒くなってから。 自分たちは手がかじかむのを必死で耐えて自転車で登校していた。 校則にはなっていないが、ポケットに手を入れて自転車に乗ってはいけないと、熱心に生活指導する先生もいる。 なぜ手袋だめなの?

ウィンドブレーカーも同様に無意味なルールで縛られていた。 黒か紺のみ、という。 いやいや、安全第一って何時も言ってるやん!とおもった。 絶対派手な色の方が自転車乗るにも、道路脇を歩くにも安全でしょ。 未だに納得いかない。

さらには、こうした校則の改正は、生徒全員が校則を守ってからじゃないとだめだ!と言い出す、お馬鹿な先生も出てきた。 どこからそんな考えが生まれてくるのだろう。 問題のある法律が存在する時、国民が全員守るまで改正すべきでない、という考え方でもあるのだろうか? とにかく生徒を平らに平らに、個性をなくし、順序付けをなくし、平等にするためなら、本末転倒もやむなし、というのが左翼先生の基本のようだった。

ソ連と社会党が大好き

授業中、露骨に自分の政治心情を生従に押しつける先生もいた。 とにかくその先生はソ連(今のロシア)と社会党(現在の社民党の母体)が大好きであった。 例えばこうである。

「ソ連は、全員が平等に扱われている。私達が見習わないといけませんね。」
「社会党も平等が理想で本当に良い政党だと、先生は思います。」
「皆さんはどう思ぃますか?」

とにかくこういう事を事あるごとに言う。 自分は、親から先生の言う事が必ずしも正しく無いので鵜呑みにしてはいけないと、釘を差されていたため全て信用した訳では無かった。 しかし、こういう聞き方をされれば小学生くらいであれば「私もそう思います」っ言えば、先生に誉められると思いそう言ってしまうものだ。 本当に信じるかどうかは別の話だが。

他にもこういう事を覚えている。 その先生が授業中にテレビをつけ、湾岸戦争のニュースを見せた。 当時は、自衛隊がPKO参加すべきか否かで議論された結果、結局参加出来ないまま、戦争が始まってしまった直後ことである。 公平に言って、生の二コースを見せて自分たちが勉強している内容との対応を示すのは非常に良い事だと思う。 しかし、先生はその後も露骨に偏った事を言う。

「自民党は、この戦争に参画しようとしてたのだから本当に悪い政党ですね。」
「PKOに反対する社会党は、戦争を阻止しようとしていますね。」
「皆さんはどう思ぃますか?」

遠足はヤマギシズム教育

遠足の行き先はヤマギシだった。正確には、「ヤマギシズム社会実顕地」(通称、ヤマギシの村)というらしい。 記憶が定かでないので、正確にこれかどうかは分からないが、バスに乗って連れて行かれ、 当地の担当者が「村」という呼び方をしていたのは覚えている。 ヤマギシをご存じない方のために簡単に説明する。ヤマギシは要するに、自由な国日本国の中にある共産主義組織である。 共産主義なので全ての人は平等に暮らすというのが理想だが、入会と同時に財産没収、 その後は「ヤマギシズム」という考え方を教育され、給料なし休日なしで農業に従事させられる。 政治団体というよりどちらかとカルト団体と見られることが多く、Amazonで「ヤマギシ会」と検索しても 「洗脳」、「カルト」などというネガティブな言葉が並ぶ。

公平に言うために、ヤマギシの良い印象も伝えておかねばなるまい。 ヤマギシが頻繁に農産物や加工品をトラックで売りに来る時期があった。 風の谷のナウシカの音楽をかけながらくるのはそのイメージが近いからだと思った。 ヤマギシが売るパンやシュークリームは当時のスーパーなどでは売ってなかったようなクオリティでめちゃくちゃ美味かったことを記憶している。今では結構あたりまとなったがシュークリームはでかくクリームが濃厚で美味かった。 天然酵母パンも当時のパン屋には売ってなかったような、手作り感あふれる味で、明らかに群を抜いて美味かった。

ヤマギシの村での体験は、前半は楽しかった。 ヤマギシでは、鶏も通常の養鶏所のようにただ餌を与えられて卵を生むだけではなく、 囲いの中を放し飼いにされていて、オスもメスも一緒くたに育てられていた。 卵はあちこちで生むので探しながらとらないといけない。 より自然な形で育てられているから良いのだという。

しかし、後半見せられたものには、今でも鮮明に記憶するほどの違和感を感じた。 1つの部屋に生徒全員が集められ、ヤマギシの村がどういうところかという説明を受けた。 まるでNHK教育テレビで出てくる体操のお兄さんのような若いお兄さんが登場して挨拶した。

「ハレハレー」

顔の両側に両手を広げて手のひらをこちら側に向け、カサゴのようにヒレをひらひらさせながら「ハレハレー」と声を放つのである。

「皆さんも一緒に!せぇ~のハレハレー」
「みなさ~~ん、ヤマギシでは喧嘩をする人は一人も居ません。誰もが笑って暮らしています。みんないつも笑顔です。」

はっきり言ってドラマ「トリック」に出てくるカルト団体そっくりである。 それ以外も色々説明してたけど、殆ど覚えていない。 この喧嘩をしないという点に子供ながら、嘘だと思ったことをはっきり覚えている。 根拠はないけど喧嘩しないやつがいる訳がない。 喧嘩の絶えない家族に生まれた自分が特殊なのかもしれないけど、とにかく率直にそう思った。

しかし、まぁ今思えば、自分の子供がこのような遠足に連れて行かれると知ったら反対するだろうな。 オオム真理教に遠足に行きますと言って怒らない親がいるだろうか? 自分の思想は自由だけど子供に偏ったことだけを見せる事自体が洗脳と同じじゃないか。 しかし、当時はそれがまかり通ってた訳だから一人の特殊な左翼先生が居ただけではないということが、今になって分かる。

反戦教育のためなら暴力もやむなし

中学時代の修学旅行は当然、左翼先生が大好きの広島になった。 反戦教育ができるからである。 確か初日は広島の手前の町で観光し、翌日大きなホールで被爆者が講演する体験談を聞いた。

中学生にとっては全く面白くないけど、それは元気盛りの中学生、友達同士が同じ宿に泊まるというだけで楽しくて仕方ない。 当然、どこの学校でも同じ、前日の夜は興奮して全く眠れない。 予め枕投げは禁止と釘を刺されたけど、たった1つのゲームを奪われただけで、中学生の自由は奪えない。 先生が居ない時は寝ないでふざけまわり、先生が来ると寝たふりをする遊びで十分楽しめた。

ほぼ寝れないまま被爆者の体験談を聞く時間になった。 薄暗いホールで座らされ、特にしゃべり上手でもない被爆者の話を聞かされれば、どんな人でも眠たくなるのと思う。 自分は真ん中からやや後ろに座っていたのだけど、前の方の男の子も女の子も多くの首が傾いているのが見えた。

自分も、うとうととし始めたその時、周りの生徒がざわめき初めて目が覚めた。 生活指導などで生徒に暴力を振るうことで有名なあの先生が寝ている生徒めがけて突進してきた(ように見えたが実際は歩いてきた)。 体も手も丸々太った190cmくらいの先生が生徒の頬を順番に バッチーーン、バッチーーン、バッチーーンと会場に音が鳴り響くくらいに連続ビンタしているのが見えた。 この先生、同じクラスのやんちゃな生徒の顔をめちゃくちゃに殴って本当に顔の形が変わるくらいボコボコにしたことがある。 実際の殴られた生徒の顔を見たが両目が腫れて戦った後のボクサー見たいになっていた。 担任の先生が大丈夫かと心配した程である。 そんな先生に連続ビンタされるのは嫌だから、その後も必死で寝ないことだけを考えて、長い講演を乗り切った。

中学生ながら後から考えてどう考えてもおかしいとおもった事を覚えている。 反戦教育は被爆者の方が可愛そうだと同情するための教育ではなくて、 暴力自体がダメなことなんだとと教えたかったに違いない。 にも関わらず、その教育をやってる目の前で教師が暴力を振るってどうするんだと。 しかも、それは警告射撃なしの即射殺のような連続ビンタである。 ヘタすれば死んでしまう生徒だっているかもしれない。 一体何を考えているのか、この先生はと思った。

こんな先生方が今義務教育組織の長になっている

これまで紹介した先生は一人ではない。 その中には、自分の得意分野を見つけてくれた、自分の人生にとって大事な大事な恩師も含まれる。 それでもやはり、偏った教育がなされていたように思う。 なぜ、義務教育の先生が特に左に偏っていたのか? 日教組の影響が強いかもしれない。 平等が正義だと考える左翼先生にとっては、とにかく他の都合の悪い事実は無視して、 というより勉強不足で知りもせず、平等であればすべてが良いと信じ込んでいた。 しかし、少し勉強すれば、事実をより多く知れば、物事はそれほど単純ではないということが分かる。 心配なのは、ここで紹介した先生方は、当時30代後半から40代前半、今は丁度教育組織の長になっている。 こうした先生が世代交代により総入れ替えするのは、50年近くかかる。 教育システムを真剣に考えるということは、本当に大事だと痛感する。

参考

2013年11月14日木曜日

Tracking関連論文のエレガントな変数定義例

PWP3D

Prisacariu, Victor and Reid, Ian: "PWP3D: Real-Time Segmentation and Tracking of 3D Objects," International Journal of Computer Vision, vol. 98, no. 3, pp. 1-20, 2012.
Let the image be denoted by $I$ and the image domain by $\Omega\subset\mathbb{R}^2$. An image pixel $x = [x,y]$ has a corresponding image value $I(x) = y$ (in our experiments an RGB value), a 3D point $X = [X,Y,Z]^T = RX_0 +T \in \mathbb{R}^3$ in the camera coordinate frame and a point $X_0 = [X_0,Y_0,Z_0]^T \in \mathbb{R}^3$ in the object coordinate frame. $R$ and $T$ are the rotation matrix and translation vector respectively (representing the unknown pose) and are parameterised by 7 parameters (4 for rotation and 3 for translation), denoted by $\lambda_i$. We use quaternions to represent rotation, so we use 4 parameters to capture the 3 degrees of freedom of the rotation.
We assume the camera calibration parameters (for each camera) to be known. Let $(f_u,f _v)$ be the focal distance expressed in horizontal and vertical pixels and $(u_o, v_o)$ the principal point of the camera. In the case of multiple cameras, the extrinsics were obtained using the Matlab Calibration Toolbox.

DTAM

Newcombe, Richard A. and Lovegrove, Steven J. and Davison, Andrew J.: "DTAM: Dense tracking and mapping in real-time", Proc. Int. Conf. on Computer Vision (ICCV2011), pp. 2320-2327, 2011.
We refer to the pose of a camera $c$ with respect to the world frame of reference $w$ as \[ T_{wc} = \left( \begin{array}{cc} R_{wc} & c_{w}\\ 0^T & 1 \end{array} \right) \] where $T_{wc} \in \mathbb{SE(3)}$ is the matrix describing point transfer between the camera's frame of reference and that of the world, such that $x_w = T_{wc}x_c$. $R_{wc} \in \mathbb{SO(3)}$ is the rotation matrix describing directional transfer, and $c_w$ is the location of the optic center of camera $c$ in the frame of reference $w$. Our camera has fixed and pre-calibrated intrinsic matrix $K$ and all images are pre-warped to remove radial distortion. We describe perspective projection of a 3D point $x_c = (x, y, z)^T$ including dehomogenisation by $\pi (x_c) = (x/z, y/z)^T$. Our dense model is composed of overlapping keyframes. Illustrated in Figure 1, a keyframe $r$ with world-camera frame transform $T_{rw}$, contains an inverse depth map $x_{r}: \Omega \rightarrow \mathbb{R}$ and RGB reference image $I_r : \Omega \rightarrow \mathbb{R}^3$ where $\Omega \subset \mathbb{R}^2$ is the image domain. For a pixel $u := (u, v)^T \in \Omega$, we can back-project an inverse depth value $d = \xi (u)$ to a 3D point $x = \pi^{-1}(u, d)$ where $\pi^{-1}(u, d) = \frac{1}{d} K^{-1}\dot u$. The dot notation is used to define the homogeneous vector $\dot u := (u, v, 1)^T$.

KinectFusion

Newcombe, Richard A. and Davison, Andrew J. and Izadi, S. and Kohli, P. and Hilliges, Otmar and Shotton, J. and Molyneaux, David and Hodges, Steve and Kim, David and Fitzgibbon, A.: "KinectFusion: Real-time dense surface mapping and tracking," Proc. 10th IEEE Int. Symp. on Mixed and Augmented Reality (ISMAR2011), pp. 127-136, 2011.
We represent the live 6DOF camera pose estimated for a frame at time $k$ by a rigid body transformation matrix: \[ T_{g,k} = \left[ \begin{array} R_{g,k} & t_{g,k}\\ 0 & 1 \end{array} \right] \in \mathbb{SE}^3 \] where the Euclidean group $\mathbb{SE}^3 := \{R, t | R \in \mathbb{SO}^3, t \in \mathbb{R}^3g\}$. This maps the camera coordinate frame at time $k$ into the global frame $g$, such that a point $p_k \in \mathbb{R}^3$ in the camera frame is transferred into the global co-ordinate frame via $p_g = T_{g,k} p_k$. We will also use a single constant camera calibration matrix $K$ that transforms points on the sensor plane into image pixels. The function $q = \pi(p)$ performs perspective projection of $p \in \mathbb{R}^3 = (x,y,z)^T$ including dehomogenisation to obtain $q \in \mathbb{R}^2 = (x/z,y/z)$. We will also use a dot notation to denote homogeneous vectors $\dot u := (u^T|1)^T$

2013年11月13日水曜日

Gnuplotによる区分関数の最小二乗法

不等号で定義域を分けて異なる関数を組み合わせた作った関数を区分関数(Piecewise Function)という。 例えば絶対値を用いた関数$f(x)=|x|$も$x=0$を境に$x>=0$における$f(x)=x$と、$x<0$における$f(x)=-x$の組み合わせと考えられる。 誤差を含むデータ点から区分関数のパラメータを求める最適化はGnuplotのfitコマンドで意外と簡単にできる。 面白いのは区分関数の区切り位置(2つの関数の境目)自体もパラメータとすることでその位置も求められること。

V字関数のフィッティング

データの作成が容易なので、次の簡単な区分関数のフィッティング方法を紹介する。 正確な名前は分からないがここではV字関数と仮名を付けておく。 この関数どういう形かというと、$a<0,\ b>0$のとき$x$軸上の$x=c$で$y=0$となり、折り返すまさにV字の関数となる。

\[ f(x)= \left\{ \begin{array}{ll} a(x-c) & (x < c)\\ b(x-c) & (x \geq c)\\ \end{array} \right. \]

上記関数に対して誤差を含むデータ点が与えられたとき、 Gnuplotのfitコマンドでパラメータ$a,b,c$を同時に非線形最適化して求めてみる。 使用したGnuplotのコマンドとデータは以下のとおり。

Gnuplotコマンド

reset

#関数の定義
f(x,a,b,c)= (x<c) ? a*(x-c) : b*(x-c)
 
#初期値の設定
a=0
b=0
c=0

#関数のフィッティング
fit f(x,a,b,c) 'sample.txt' via a,b,c
 
#表示設定
set xr[-10:10]
set samples 100
set terminal wxt

#関数とデータの描画
plot f(x,a,b,c) 
replot 'sample.txt' lc 3

データ "sample.txt"

-4 11
-3 7.5
-2 6
-1 3
-0 2
1 0
2 5
3 12
4 16
5 20
6 27

このデータは、$a=-2, b=5, c=1$の時の関数の値に対してわずかな誤差を人手により与えている。 最適化の初期値は意図的に正解値から遠ざけて$a=0, b=0, c=0$としたが、おおよそ上の正解値に近づいた結果$a=-2.0544, b=5.21874, c=0.936731$が得られている。

データとフィッティングされた関数

フィッティング結果

gnuplot> load 'fitting_sample.txt'
 
 Iteration 0
 WSSR        : 1780.25           delta(WSSR)/WSSR   : 0
 delta(WSSR) : 0                 limit for stopping : 1e-005
 lambda   : 1.91485
 
initial set of free parameter values
 
a               = 1e-030
b               = 1e-030
c               = 1e-030
/
   :
  (省略)
   :
 
 Iteration 6
 WSSR        : 5.75979           delta(WSSR)/WSSR   : -9.82704e-009
 delta(WSSR) : -5.66017e-008     limit for stopping : 1e-005
 lambda   : 1.91485e-005
 
resultant parameter values
 
a               = -2.0544
b               = 5.21874
c               = 0.936731
 
After 6 iterations the fit converged.
final sum of squares of residuals : 5.75979
rel. change during last iteration : -9.82704e-009
 
degrees of freedom    (FIT_NDF)                        : 8
rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 0.848513
variance of residuals (reduced chisquare) = WSSR/ndf   : 0.719974
 
Final set of parameters            Asymptotic Standard Error
=======================            ==========================
 
a               = -2.0544          +/- 0.1336       (6.501%)
b               = 5.21874          +/- 0.1976       (3.787%)
c               = 0.936731         +/- 0.1152       (12.3%)
 
correlation matrix of the fit parameters:
 
               a      b      c      
a               1.000 
b               0.403  1.000 
c               0.490  0.822  1.000・

参考

2013年11月9日土曜日

cv::gpu::flipの入出力に同じ変数は使えない

OpenCVのGpuMatは非常に便利だけど、GPGPUの仕組みを全く理解せずにBlack Boxとして使うと時々問題になる。 例えばcv::gpu::flip()関数。 CPUのみを使う関数cv::flipでは入出力に同じ変数を用いても問題にならなかったが、 cv::gpu::flip()の場合は正常にflipされない(OpenCV2.6.4で確認)。 次の有名な標準画像を入力として試してみた。

#include <opencv2/opencv.hpp>
#include <opencv2/gpu/gpu.hpp>
 
void main(int argc, char *argv[])
{
    cv::Mat image = cv::imread("lena_std.png");
    cv::gpu::GpuMat image_gpu;

    //画像をCPUからGPUに送る
    image_gpu.upload(image);
 
    //CPU上で画像の上下反転、保存
    cv::flip(image, image, 0);
    cv::imwrite("fliped_lena_cv.png", image);
 
    //GPU上で画像の上下反転(srcとdstを同じにしている)
    cv::gpu::flip(image_gpu, image_gpu, 0);
 
    //画像をGPUからCPU側に戻しファイルに保存
    image_gpu.download(image);
    cv::imwrite("image_gpu.png", image);
}

上記コードで出力される画像は次のとおり。 左はcv::flip()、右はcv::gpu::flip()の結果。 GPUの方は、正常に画像が生成されていないことが分かる。

次のコードは上記コードを少し修正して、GpuMatを2つ用意して、 cv::gpu::flip()の入出力の変数を変えてみた。

#include <opencv2/opencv.hpp>
#include <opencv2/gpu/gpu.hpp>

void main(int argc, char *argv[])
{
    cv::Mat image = cv::imread("lena_std.png");
    cv::gpu::GpuMat image_gpu;
    cv::gpu::GpuMat image_gpu2;

    //画像をCPUからGPUに送る
    image_gpu.upload(image);
 
    //CPU上で画像の上下反転、保存
    cv::flip(image, image, 0);
    cv::imwrite("fliped_lena_cv.png", image);
 
    //GPU上で画像の上下反転(srcとdstを変えた)
    cv::gpu::flip(image_gpu, image_gpu2, 0);
 
    //画像をGPUからCPU側に戻しファイルに保存
    image_gpu2.download(image);
    cv::imwrite("image_gpu22.png", image);
}

左はcv::flip()、右はcv::gpu::flip()の結果。 今回はCPU、GPU共に同じ結果が得られた。

2013年11月7日木曜日

OpenGLで超簡単なオフスクリーンレンダリング

やりたい事

やりたい事はタイトルの通り。 GLUTを使うとOpenGLを簡単に扱える反面、glutMainLoop()に全部制御を持って行かれてしまう。 FBO(frame buffer object)のような難しいものは触りたくない。 ウィンドウを表示せずにOpenGLで描いた画像を取得する方法を探しても意外と見つからなかったのでここにメモしておく。 要は、glutCreateWindow()関数を呼んでウィンドウを作る準備をしておけば、 バッファは利用可能な状態になっているようなので、 実際にウィンドウが表示されなくても描画したり、描画したものを画像として取得したりできるということ。 この方法だと複数のウィンドウも可能だけど、1つのウィンドウを表示してもう一つのウィンドウは表示しないということは出来ないかも。

サンプルコード

#include <gl/glut.h>
#include <opencv2/opencv.hpp>
 
void main(int argc, char *argv[])
{
    int w = 256;
    int h = 128;
    cv::Mat_<uchar> image(h, w);
    glutInit(&argc, argv);
    glutInitWindowSize(w, h);
    glutCreateWindow(argv[0]);
 
    //描画
    glClearColor(0.1, 0.1, 0.1, 1.0);
    glClear(GL_COLOR_BUFFER_BIT);
    glutSolidCube(0.5);
    glFlush();
 
    glReadPixels(0, 0, w, h, GL_RED, GL_UNSIGNED_BYTE, image.data);
 
    cv::imwrite("test.png", image);
}

出力画像

2013年10月22日火曜日

デュアルカメラ機能付きのスマホ、タブレット比較

従来スマホやタブレットの前面(内側)と後面(外側)に装備されたビデオカメラを同時に起動することができなかったが、 2013年頃からようやくその機能をもつタブレット等が売り出されるようになってきた。 この機能の名称は以下のように異なる名前で呼ばれるが、共通点はデュアル(dual)のようだ。
  • デュアルカメラ(dual camera)
  • デュアルショット(dual shot)
  • デュアル録画(dual video recording)
このページではデュアルカメラに統一し、2013年10月21日現在この機能をもつタブレット、スマホをまとめた。

GALAXY Note3 SC-01F

Samsung GALAXY S4-04E

LG Optimus G pro L-04E

参考

2013年5月11日土曜日

スターウォーズを見れば分かる敵国日本

スターウォーズ大好き

スターウォーズは、大好きな映画の1つ。 特に、古い3作、エピソード4、5、6はハリウッド映画の最高傑作だと思う。 子供の頃テレビで放送されたものをVHSに録画して擦り切れるほど何度も見た。 宇宙空間をものすごいスピード飛び回る宇宙船は、全人類の男の子共通の夢じゃないか。 未だに、ミレニアムファルコン号が狭い隙間を縫って飛び回るシーンを想像することがある。

スターウォーズの帝国軍=大日本帝国軍

そんな大好きなスターウォーズだけど、最近この映画が第二次世界大戦、特に日本そのものを描いていることに気づいた。 第二次世界大戦中の日本と言えばもちろん大日本帝国。 以下、スターウォーズと大日本帝国との共通点をあげてみる。

皇帝天皇
ダースベーダーの姿甲冑姿
ライトセイバー日本刀
フォース精神力
同盟軍連合軍
スーパー・スター・デストロイヤー戦艦大和
突進してくる帝国軍の兵隊特攻してくる日本兵

どうだろう、この偶然の一致(笑)。 少し無理やりかもしれないけど、スーパー・スター・デストロイヤーが共和国軍の小型宇宙船に集中砲火を浴びて沈み行くシーンと、 (映画か何かで見たのだろうか)戦艦大和がやられるシーンとの類似性を発見した時には鳥肌が立った。 ストーリーはそのままに言葉とキャラクターを第二次世界大戦のそれに置き換えれば、案外普通のアメリカ視点のドキュメンタリー映画になってしまうんじゃないか、とすら思う。

染み付いた「帝国=悪」のイメージ

でも本当に偶然なのか? 「フォース」や「ジェダイ」などの言葉は、日本の精神性を連想させるし、彼らの戦いはチャンバラそのものだ。 ジョージ・ルーカスが黒澤監督の映画に影響を受けたという話は、嘘か本当か分からないが有名な話。 これが本当だとすれば、日本人である黒澤監督を尊敬していたということであり、 日本を敵視してこのようなストーリーにしたとは考えにくい。

私の想像は、第二次世界大戦という出来事は、アメリカなど戦勝国の教育に深く心に刻み込まれ、 悪の集団の極みを想像するとすれば、それがどうしても当時の日本とナチスになってしまう、んじゃないかかということ。 たぶん日本がどんな悪いことをしたのか詳しく説明できる人は少ない(せいぜいパールハーバーの不意打ちを挙げるくらいだろう)。 それでも悪と言えば、例え日本に敬意を払う人ですら、絶対悪としてのナチスと日本が連想して出てくるんじゃないか。 そういえば、スターウォーズの帝国軍兵士の制服はナチスの制服を連想する。

「帝国=悪」のイメージは結構一般的になっているように思う。 そもそもなぜそうなんだろう? このスターウォーズのせい?天皇=大日本帝国軍の戦犯のイメージから? 「帝国主義」という言葉もある。これは完全に悪のイメージ。 それともGHQのせいで「大日本帝国=悪」と教育されたから? 「帝国ホテル」や「帝国大学」と聞いて悪のイメージはない。 帝国≒王国なのだけど「王国」と聞いて悪のイメージはない。 フランス革命以後、王様を倒して自由を勝ち取り、民主主義に変わったという西洋のイメージから、 「帝国=悪」になったのではないか?

いずれにせよ、絶対悪としてのナチスと大日本帝国を連想してくれなければ、 アメリカが非戦闘要員の住宅地に原爆を落としたり、 焼夷弾で町ごと焼き払ったことの正当性は失われてしまうのは事実だと思う。 だからスターウォーズの正義と悪の対立構造を見て、大日本帝国軍のような帝国軍を見て、 すんなりとどちらが正義でどちらが悪かを、認識できるのだと。

国連は未だに日本を敵国としている

では、大日本帝国=悪のイメージは、頭に染み付いたイメージなのか? 実際問題として、国際連合憲章に「敵国条項(Enemy Clauses)」というのがあるらしい。 ここには、2013年5月11日現在、未だに第二次世界大戦当時の敵国(日本、ドイツなど)を悪として扱う旨の記述がある。

第53条

1. 安全保障理事会は、その権威の下における強制行動のために、適当な場合には、前記の地域的取極または地域的機関を利用する。但し、いかなる強制行動も、安全保障理事会の許可がなければ、地域的取極に基いて又は地域的機関によってとられてはならない。もっとも、本条2に定める敵国のいずれかに対する措置で、第107条に従って規定されるもの又はこの敵国における侵略政策の再現に備える地域的取極において規定されるものは、関係政府の要請に基いてこの機構がこの敵国による新たな侵略を防止する責任を負うときまで例外とする。

2. 本条1で用いる敵国という語は、第二次世界戦争中にこの憲章のいずれかの署名国の敵国であった国に適用される。

直訳だからだろうか、意味不明の文章だけど、よく読むと要するにこう書かれてる。 まず1で、軍事力を使うなら国連の許可をとれよと。もっとも、敵国は例外だけど。 次に2で、敵国ってようするに第二次世界大戦中の敵国のこと。 つまり、この2つを組み合わせると日本への軍事行動は国連の許可なしでもOKと書かれている! 隣国などがこれを悪用しないかと心配になるのは当然でしょ。

そもそも国連という名前がおかしい

国連は国際連合の短縮形として学校で習ったが、そもそもこの訳がおかしい。 英語で書くと「United Nations」である。 直訳すれば「連合国」であり、それは第二次世界大戦で勝利した側である「連合国軍」の連合国である。 連合国軍は英語はAlliesが使われているようだが、連合国軍の目的や敵国が何処なのかを明記した「連合国共同宣言(The Declaration by United Nations)」は、日本の敗戦前に作られた文章であり、やはり連合国軍=The United Nationsであることが分かる。 決して「連合国」を「国際連合」という新しい単語にすり替える必要はない。

ズボンのことをボトムスなどと呼び変えて、新しいものかのように装うファッション業界とは違うので、 本来名前なんかどうでもいいのだが、ここに歴史的経緯がある以上、その名前がついた経緯が重要である。 誰でも知っているとおり、安全保障理事会の常任理事国は戦勝国そのものである。 そして無条件に拒否権というどう考えても差別的な権利を持っている。 その国が世界の安全保障を考えてようが、考えてまいが、である。 要するに最終的には戦勝国の都合が正義としてまかり通る組織だということだ。

民主主義を理念とする新しい共同体を

この組織の常任理事国が自分たちの特権を手放すことは、決して、無いだろう。 この不平等があるかぎり、必ずその不満がたまることになる。 常任理事国よりも国力のある国が、今は日本やドイツくらいしかないかもしれないが、 これから増えてきた時にはその不満が爆発する危険性が常にある。 最大の武器である民主主義という理念、国民一人ひとりの自由を大切にするという原則を持った国どおしが結束し、 平等な権利でその次代に国力のある国がリーダーシップを取れるような仕組みをもつ、 「国際民主主義連合」のようなものを作っていくべきではないかと思う。

スターウォーズのちょっとした発見から、ここまで連想するのはやり過ぎだろうか?

参考

2013年3月30日土曜日

Sony Vaio Duo 11 Review: AR developers' perspective

This review is only interested in whether the Sony Vaio Duo 11, a Window 8 tablet PC, can be used as an augmented reality (AR) terminal. The most important advantages of Windows 8 tablet as far as this review are that windows developers are not forced to change their programming environment, and that Windows 8 devices only offer the option of using Intel CPUs. Therefore, this review does not mention the usability of either normal Windows 8 GUI or pre-installed software.

Configuration related to AR

  • OS: Windows 8 Pro 64
  • CPU: Intel Core i7-3687U (max 3.30GHz, 2 cores, 4 threads)
  • GPU: Intel HD Graphics 4000
  • RAM: 8GB
  • rear camera: Full HD web camera
  • front camera: Full HD web camera
  • pose sensor:accelerometer,gyro,magnetic field sensor
  • weight: 1.665 kg
  • size: 319.9mm×199mm×17.85mm

Limitation on Camera Capture

If you use Vaio Duo 11 as an AR device, the camera, together with the display, becomes an indispensable part. The Vaio duo 11 has front and rear cameras which are placed in opposite directions to each other. The rear camera can be used to obtain a real scene to be displayed with virtual (CG) objects. The front camera can also be used to detect user's action for interaction in augmented reality space. In this sense, it is interesting whether or not both cameras can work simultaneously, which has not yet been seen. The following list shows the possible combination of video resolutions.

  • front: 320x240 (7.5fps), rear: 1280x1024 (7.5fps) (frame droping)
  • front: 640x480 (30fps), rear: 640x480 (30fps) (frame droping)
  • front: 640x480 (30fps), rear: 320x240 (30fps)
  • front: 320x240 (30fps), rear: 640x480 (30fps)
  • front: 320x240 (30fps), rear: 320x240 (30fps)

Possible choice in practical use for AR is 640x480 pixels for one camera and 320x240 for the other. What is surprising is that each camera cannot capture full HD resolution video (1920 x 1080) at 30 fps. "Full HD web camera" means full HD sill camera. This is a really misleading expression. The above list were obtained by testing the following OpenCV code. DirectShow GraphEdit did not work well.

#include <opencv2/opencv.hpp>
#include <iostream>
using namespace cv;
 
int main(int, char**)
{
    VideoCapture cap1(0);
    cap1.set(CV_CAP_PROP_FRAME_WIDTH, 320.0);
    cap1.set(CV_CAP_PROP_FRAME_HEIGHT,240.0);
 
    VideoCapture cap2(1);
    cap2.set(CV_CAP_PROP_FRAME_WIDTH,640);
    cap2.set(CV_CAP_PROP_FRAME_HEIGHT,480);
 
    std::cout &lt;&lt; cap1.get(CV_CAP_PROP_FPS);
    std::cout &lt;&lt; cap2.get(CV_CAP_PROP_FPS);
 
    cap1.set(CV_CAP_PROP_FPS,30);
    cap2.set(CV_CAP_PROP_FPS,30);
 
    double sum_t = 0;
    for(int i=0;;i++)
    {
        double t = (double)getTickCount();
        Mat frame1;
        cap1 &gt;&gt; frame1;
        imshow(&quot;cap1&quot;, frame1);
 
        Mat frame2;
        cap2 &gt;&gt; frame2;
        imshow(&quot;cap2&quot;, frame2);
 
        sum_t += getTickFrequency()/((double)getTickCount() - t);
        if(i % 10 ==0)
        {
            std::cout &lt;&lt; sum_t/10 &lt;&lt; &quot;\r&quot;;
            sum_t =0;
        }
        if(waitKey(1) &gt;= 0) break;
    }
    return 0;
}

2013年3月7日木曜日

PTAM/PTAMMの変数アクセス法

Trackerクラスからのアクセス

現在のキーフレームで観測される特徴点の三次元位置を列挙したい

PTAM/PTAMMの特徴はキーフレームという代表フレーム上の特徴点のみが三次元復元される。 現在のフレームに最も違いカメラ位置であるキーフレームで観測される特徴点のみを列挙したい場合には、 Trackerクラスのメンバ変数mCurrentKFから特徴点を参照すれば良い。

KeyFrame & nearestKF = mMapMaker.ClosestKeyFrame(mCurrentKF); //もともとprivateだがpublicに変更した
std::map::iterator p;
std::map::iterator s = nearestKF.mMeasurements.begin();
std::map::iterator e = nearestKF.mMeasurements.end();
for(p = s; p != e; ++p)
{
 MapPoint& mp = *(p->first);
 if(mp.bBad) continue; //三次元位置が推定されていない点も登録されている
 std::cout 
  << mp.veWorldPos[0] << "," 
  << mp.veWorldPos[1] << "," 
  << mp.veWorldPos[2] << std::endl;
}

現フレームでトラックされた特徴点の三次元位置を列挙したい

現フレームで新たに検出された特徴点のうち、 既に三次元位置が登録されている特徴点とのマッチングを終え、 マップ点の三次元位置と関連付けられている特徴点のみを列挙したい場合は、 特徴点を参照するタイミングに注意せねばならない。 単純には、以下のmy_process()ようにトラッキングを終えた直後に、参照すればよい。

void Tracker::TrackFrame(Image<CVD::byte> &imFrame, bool bDraw)
{
    :
    :
 
    // Decide what to do - if there is a map, try to track the map ...
    if(mpMap->IsGood())
    {
        if(mnLostFrames < NUM_LOST_FRAMES)  // .. but only if we're not lost!
        {
            if(mbUseSBIInit)
                CalcSBIRotation();
            ApplyMotionModel();       // 
            TrackMap();               //  These three lines do the main tracking work.
            UpdateMotionModel();      // 
 
            //このタイミングでトラッキングされたばかりの特徴点データを取得する。
            my_process();

この箇所では、TrackerDataからとランキングされている特徴点の三次元位置を以下のように取得する。

void Tracker::my_process()
{
    vector<TrackerData*>::iterator it = vIterationSet.begin();
    for(; it!= vIterationSet.end(); ++it)
    {
        if(! (*it)->bFound) continue;
 
        //2d
        Vector<2> & pos2d = (*it)->v2Found;
        cout << boost::fromat("2d point: %f, %f\n") % pos2d[0] % pos2d[1];
 
        //3d
        Vector<3> & pos3d = (*it)->Point.v3WorldPos;
        cout << boost::fromat("3d point: %f, %f, %f\n") % pos3d[0] % pos2d[1] % pos2d[2];
    }
}

2013年1月30日水曜日

拡張現実用 Windows8 高性能タブレット 比較

高度な計算ができる拡張現実感(AR)用途の端末に必要な機能に焦点を絞ってWindows 8タブレットを比較してみた。 CPU Intel Core i7は必須とし、オプションは可能な限りハイスペックなものを選ぶ。 (2013/1/30)

  • OS
  • CPU
  • GPU
  • メモリ
  • 外側カメラ
  • 内側カメラ
  • 姿勢センサ
  • 重さ
  • サイズ(幅×高さ×厚さ)

Sony VAIO Duo 11

  • OS: Windows 8 Pro 64
  • CPU: Intel Core i7-3687U (max 3.30GHz, コア2, スレッド4)
  • GPU: Intel HD Graphics 4000
  • メモリ: 8GB
  • 外側カメラ: Full HDウェブカメラ
  • 内側カメラ: Full HDウェブカメラ
  • 姿勢センサ:加速度センサー,ジャイロ,地磁気センサー
  • 重さ: 1.665 kg
  • サイズ: 319.9mm×199mm×17.85mm

DELL XPS 12 Ultrabook

  • OS: Windows 8 Pro 64
  • CPU: Intel Core i7-3537U (max 3.10GHz, コア2, スレッド4)
  • GPU: Intel HD Graphics 4000
  • メモリ: 8GB
  • 外側カメラ: -
  • 内側カメラ: 1.3 MP Webカメラ、1,280 x 1,024ピクセル
  • 重さ: 1.52kg
  • サイズ: 317mm×215mm×20mm

Panasonic Let's note AX2

  • OS: Windows 8 Pro 64
  • CPU: Intel Core i7-3517U (max 3.00 GHz, コア2, スレッド4)
  • GPU: Intel HD Graphics 4000
  • メモリ: 4GB
  • 外側カメラ: -
  • 内側カメラ: 720Pウェブカメラ
  • 姿勢センサ:地磁気センサー, ジャイロセンサー, 加速度センサー,照度センサー
  • 重さ: 1.14kg
  • サイズ: 288mm×194mm×19mm

Lenovo IdeaPad Yoga 13

  • OS: Windows 8 64bit
  • CPU: Intel Core i7-3517U (max 3.0GHz, コア2, スレッド4)
  • GPU: Intel HD Graphics 4000
  • メモリ: 8GB
  • 外側カメラ: -
  • 内側カメラ: HD 720p
  • 重さ: 1.5kg
  • サイズ: 333.4mm×224.8mm×16.9mm

2013年1月17日木曜日

CSDPのVisual Studio2010でのコンパイル方法

CSDPは半正定値計画問題(Semidefinite Programming)を取り扱うC言語ライブラリ。 内部ではLPACKが使用されており予めインストールされている必要がある。 ここでは、CLAPACKをソースコードからコンパイルして準備する方法も含めて書き留めておく。

動作環境

  • Visual Studio 2010 SP1
  • Windows 7 Pro SP1

CLAPACKをインストール

/MTオプションでコンパイルされたものや、Intel Compilersが必要なものしか見つからなかったので、ソースコードからコンパイルすることにした。Visual Studio用のパックが準備されているが、そのままではコンパイルが通らない。プロジェクトの設定を変更してもすべてのエラーを取り除くことはできないが、ここではCSDPに使うことのみを考え、blas, clapack, libf2cの3つのプロジェクトのみコンパイルする。

  1. CLAPACK-3.1.1-VisualStudio.zipをダウンロード。
    http://www.netlib.org/clapack/
  2. clapack.slnを開いてVisual Studio 2010形式に変換する。
  3. BLAS/SRC内の以下のファイルをblas.vcxprojに追加する。
    • csrot.c
    • drotm.c
    • drotmg.c
    • dsdot.c
    • sdsdot.c
    • srotm.c
    • srotmg.c
    • zdrot.c
  4. ビルドモードをRelease without wrapにする。
  5. blas, clapack, libf2cのプロジェクトプロパティ/コード生成/ランタイムライブラリでマルチスレッド(/MT)をマルチスレッドDLL(/MD)に変更する。
  6. libf2cプロジェクトのfio.hの111行目をコメントアウトする。
    //extern int isatty(int);
    
  7. blas, clapack, libf2cプロジェクトのみビルド。

CSDPをインストール

  1. CSDPのページにアクセスし、 「​compressed tar archive」からCsdp-6.1.1.tgzをダウンロードして解凍。
  2. ここから Csdp-6.1.1_vcproj.zipをダウンロード。
  3. 解凍するとlibsとvcbuildというフォルダができるのでこの二つをCsdp-6.1.1フォルダ直下にコピー。
  4. vcbuild内のvcbuild.slnを開いてビルドする。

参考

2013年1月12日土曜日

PWP3Dの動かし方

画像から三次元物体の位置姿勢を6自由度で推定するPWP3Dのソースコードが公開されている。1つの物体用と複数の物体用の2種類のコードが公開されているが、そのうちPerseusLib_beta2のRelease Win32をコンパイルして動作確認する。これもgSLICと同様にプロジェクトファイルの問題がありすぐには動作しなかったので対処法を書き残す。ここではCUDAはインストール済みとする。

動作環境

  • Visual Studio 2010 SP1
  • Windows 7 Pro SP1
  • Nvidia Quadro 4000 for Mac
  • NVIDIA CUDA Toolkit v4.2(インストール済み)

PWP3Dのコンパイル

  1. ここからPWP3D [beta 2 (multiple objects, multiple (independent) views) - source - VS2010 64bit](perseuslib_beta2.zip)をダウンロードする。
  2. PerseusLib.slnをVisual Studio 2010に変換する。
  3. gSLICの動かし方」と同様にPerseusLib.vcprojの設定の問題でコンバートできないことを確認する。
  4. PerseusLib.vcxprojの45行目と257行目、Perseus.vcxprojの45行目と195行目にある「CUDA 4.0.props」を「CUDA 4.2.props」に修正する。
    修正前
      <ImportGroup Label="ExtensionSettings">
        <Import Project="$(VCTargetsPath)\BuildCustomizations\CUDA 4.0.props" />
      </ImportGroup>
    
    修正後
      <ImportGroup Label="ExtensionSettings">
        <Import Project="$(VCTargetsPath)\BuildCustomizations\CUDA 4.2.props" />
      </ImportGroup>
    
  5. Perseus.slnを起動するとPerseusとPerseusLibプロジェクトが正しくロードされるので、Perseusプロジェクトの追加のインクルードディレクトリに書かれているライブラリの場所を消す。これは消さなくても問題ないが、全く必要なくややこしいので消すことにする。
    C:\CUDA\v3.2\include;
    C:\SDK\FreeImage\include;
    C:\SDK\Qt\4.7.4\include;
    C:\SDK\CLEyeMulticam\include;
    C:\SDK\VideoInput;
    
  6. FreeImage version 3154 Win32をインストール。ライブラリを正しく設定。
  7. PerseusLibプロジェクトのOther/Header Files/PerseusLibDefines.hの3~7行目をコメントアウト。
    #ifdef NDEBUG
    #pragma comment( lib, "C:\\SDK\\FreeImage\\x64\\lib\\FreeImage.lib" )
    #else
    #pragma comment( lib, "C:\\SDK\\FreeImage\\x64\\lib\\FreeImage.lib" )
    #endif
    
  8. PerseusLibプロジェクトの CUDA/Header Files/CUDADefines.hの3、4行目をコメントアウト。
    #pragma comment( lib, "C:\\CUDA\\v4.0\\lib\\x64\\cuda.lib" )
    #pragma comment( lib, "C:\\CUDA\\v4.0\\lib\\x64\\cudart.lib" )
    
  9. PerseusLibプロジェクトのプロパティからCUDA C/C++ / Device / Code Generationの以下の値を変える。
    修正前
    compute_10,sm_10
    
    修正後
    compute_20,sm_20
    
  10. Perseusプロジェクトの追加のライブラリにを追加。
    cudart.lib;cuda.lib;FreeImage.lib;
    
  11. PerseusおよびPerseusLibプロジェクトのプリプロセッサに_ITERATOR_DEBUG_LEVEL=0を追加。

参考

2013年1月4日金曜日

gSLICの動かし方

リアルタイムで画像を小さな領域(superpixel)に領域分割する手法gSLICのコードをテストしてみた。ソリューションファイル(*.sln)をロードする段階でプロジェクトファイル(*.vcxproj)が正常に読み込まれなかった。CUDAのインストールやプロジェクトの設定で追加しなければならなかったオプションなどを記録しておく。(2013/01/04)

動作環境

  • Visual Studio 2010 SP1
  • Windows 7 Pro SP1
  • Nvidia Quadro 4000 for Mac
  • NVIDIA CUDA Toolkit v4.2

CUDA4.2をインストール

CUDAを利用したgSLICのプログラムをコンパイルするためには、CUDA Toolkit、CUDA Driver、GPU Computing SDKが必要なのでこれらを順にインストールする。

  1. CUDA 4.2 Toolkit 32bitをダウンロード。64bitマシンでも32bitコンパイルするなら32bitをダウンロードする。
    https://developer.nvidia.com/cuda-toolkit-42-archive
  2. ダウンロードされたcudatoolkit_4.2.9_win_32.msiを実行してインストール。Typicalを選択する以外特殊な操作はない。
  3. CUDA 4.2 Driver 64bit()をダウンロード&インストール。32bitコンパイルする場合でも、64bitマシンならドライバ64bitを選ぶ。
    https://developer.nvidia.com/cuda-toolkit-42-archive
  4. ダウンロードされたdevdriver_4.2_winvista-win7_32_301.32_general.exeを実行してファイルを解凍。 「NVIDIA\DisplayDriver\275.89\WinVista_Win7_64\International」に解凍されたsetup.exeを実行してインストール開始。既存のドライバの設定を残すオプション「高速」を選択。
  5. GPU Computing SDK4.2.9をダウンロード。
    https://developer.nvidia.com/cuda-toolkit-42-archive
  6. ダウンロードされたgpucomputingsdk_4.2.9_win_32.exeを実行してインストール。

ルールファイルをコピー

  1. C:\Program Files (x86)\NVIDIA GPU Computing Toolkit\CUDA\v4.2\extras\visual_studio_integration\rulesを開く。
  2. C:\Program Files (x86)\Microsoft Visual Studio 10.0\VC\VCProjectDefaultsも開く。
  3. NvCudaDriverApi.v4.2.rulesとNvCudaRuntimeApi.v4.2.rulesを最初のフォルダから次のフォルダにコピーする。

gSLICのコンパイル

何よりもはじめにgSLIC自身をダウンロードして来なければならないが、2013/01/04現在GoogleでgSLICと検索すると古いページがトップに出てくる。こちらはVisual Studio 2008のバージョンだがこれは古いので扱わない。新しいページからダウンロードしてVisual Studio 2010のバージョンを用いる。gSLICのデモプログラムにはCUDA以外にOpenCVも使用されている。予めOpenCVが正しくインストールされているとする。これ以外に必要な設定は以下のとおり。

  1. 新しいページからgSLIC_vs2010.zipをダウンロードする。
  2. Visual Studioを起動するとエラーがでることを確認する。
    1. ダウンロード・解凍して生成されるフォルダ内のgSLIC.slnからVisual Studioを起動すると「Microsoft Visual Studio: ソリューション内の1つ以上のプロジェクトが正しく読み込まれていません。詳細については、出力ウィンドウを確認してください。」というエラーウィンドウが表示される。
    2. 出力ウィンドウに「gSLIC\gSLIC\gSLIC.vcxproj : error : インポートされたプロジェクト "C:\Program Files (x86)\MSBuild\Microsoft.Cpp\v4.0\BuildCustomizations\CUDA 4.0.props" が見つかりませんでした。 宣言のパスが正しいかどうか、およびファイルがディスクに存在しているかどうかを確認してください。 gSLIC\gSLIC\gSLIC.vcxproj」の類のエラーが出てることを確認する。
  3. gSLIC.vcxprojファイルのCUDA関連項目を直接テキストエディタで修正する。修正箇所は4箇所。(絶対パス出ない方がよいがとりあえず・・・)
    1. 25行目。CUDAのバージョンとNVIDIA GPU Computingインストール先のパスを修正。
      修正前
      <CudaToolkitCustomDir>
          ..\..\..\..\..\..\..\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v4.0
      </CudaToolkitCustomDir>
      
      修正後
      <CudaToolkitCustomDir>$(CUDA_INC_PATH)</CudaToolkitCustomDir>
      
    2. 53行目。CUDAのバージョン修正。
      修正前
      <Import Project="$(VCTargetsPath)\BuildCustomizations\CUDA 4.0.props" />
      
      修正後
      <Import Project="$(VCTargetsPath)\BuildCustomizations\CUDA 4.2.props" />
      
    3. 90行目。同じ。
      修正前
      <AdditionalLibraryDirectories>
          C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v4.0\lib\Win32;
          %(AdditionalLibraryDirectories)
      </AdditionalLibraryDirectories>
      
      修正後
      <AdditionalLibraryDirectories>
          $(CUDA_LIB_PATH);%(AdditionalLibraryDirectories)
      </AdditionalLibraryDirectories>
      
    4. 172行目。同じ。
      修正前
      <Import Project="$(VCTargetsPath)\BuildCustomizations\CUDA 4.0.targets" />
      修正後
      <Import Project="$(VCTargetsPath)\BuildCustomizations\CUDA 4.2.targets" />
      
  4. OpenCV関連の修正。ビルドするとCUDAのビルドは正常終了するがOpenCVのライブラリの場所に問題がある旨のエラーがでるかもしれない。OpenCVのライブラリが絶対パスなどで設定されているのが原因。
    1. 追加のインクルードディレクトリに「C:\OpenCV2.1\include\opencv」が設定されているので、これを削除または正しいディレクトリを設定する。
    2. 上と同様に追加のライブラリディレクトリに設定されたC:\OpenCV2.1\libを削除する。(これは消さなくても問題ない)
    3. main.cppのインクルードを正しく修正。
      修正前
      #include "cv.h"
      #include "highgui.h"
      
      修正後
      #include <opencv/cv.h>
      #include <opencv/highgui.h>
      
    4. リンカー/入力/追加の依存ファイルにcxcore210d.lib;cv210d.lib;highgui210d.libが設定されているのでこれを削除。
    5. OpenCV ver.2 超々々入門(C++)」の方法でOpenCVのlibファイルを指定する。

参考