2014年11月18日火曜日

Juliaで外部のプログラムを動かすときのメモ

Juliaで外部のプログラムを動かすときのメモ。
PythonではPopenとかを使ってたので、Juliaでもそれと同じようなものがあるかなーと思ったけど、なかなか見つからなかった。

公式マニュアルはこれ。
http://julia.readthedocs.org/en/latest/manual/running-external-programs/

外部プログラムを起動して、そこから標準入力を扱いたいときは、こっちが参考になった。http://blog.leahhanson.us/running-shell-commands-from-julia.html

概要は、
プログラムが起動して、
write()で書き込み。書き込んだらその出力がsoにたまっていくので、
readall()でその出力を得る。
(so,si,pr) = readandwrite(`yourprogram`)
write(si, "command1")
write(si, "command2")
close(si)
output = readall(so)


これも一応メモ。
http://julialang.org/blog/2013/04/put-this-in-your-pipe/


2014年11月11日火曜日

Novel Class

アメリカ滞在の後半はNovel Classを取った。文字通り小説の授業だ。

半期に3つの小説を読んだ。1つはベトナム戦争についての話。もうひとつはThe Catcher In The Rye。最後にSidartha。Sidarthaはヘルマン・ヘッセの小説。どれも一字一句辞書で調べながらじっくり読んだんで、ものすごく思い出に残っている。

ベトナム戦争についての話はあんまり覚えてない。ただ、そのなかの参考資料にBreaking Awayという映画があった気がする。本当にこの小説の授業で出会ったものかは定かではないえけど、ものすごく感傷にひたったのを覚えている。Breaking Awayは青春を舞台にしたもので4人の登場人物がそれぞれ精一杯行きている。あーおれはこんなところで独りでなにをやっているんだろう。と心の通った友達とのチームプレイがいかに大事かが身に染みたのを覚えている。

2つ目のThe Catcher In the Ryeは有名な小説だが、村上春樹の日本語訳はまったくおすすめしない。帰国後日本語訳を手にとって読んでみたけど、まったく感情移入することが出来んかった。こっちはあまり前向きな話ではない。最後にDon't tell anybody anything. You will miss somebody.というのがすごく思い出に残っている。他にもI have to admit itとかThat kills me.とかものすごく当時のこころに響くセリフばかりだった。

Sidarthaは難しい本が読めなかったので、優しい英訳の本を読んだ。Sidarthaが悟りをひらくまでの苦悩が思い出に残っているが、詳しい内容は忘れたw


それぞれの本を読んだ後の課題は自由。どんな方法でも読んだ感想を表現すればよい。模型をつくった人もいれば、冊子を作った人もいた。日本のいわゆる読書感想文みたいに決まった書式はない。自分の得意な方法で表現すればよかった。いま日本にこんな教育をするところはあるんかな?


2014年11月6日木曜日

Julia Array操作の基本をメモ。

Pythonでやっていたlistの操作にすっかり慣れてしまったけど、JuliaのArrayではどうやるか確かめる。
Julia0.3で試して動いたものを掲載。

配列の生成

a = [i for i in 1:10]     #1から10を要素に持つ配列生成
a = [1:10]                #同上
b = reshape([1:9], (3,3)) #3x3の2次元配列。
b = [1 2 3; 4 5 6; 7 8 9] #3x3の2次元配列。上とは行、列が逆(転置の関係)

要素数

length(a)

要素アクセス

a[1]
a[1:end]   #要素1から最後まで
a[2:3:8] #要素2から8を3個おきに

追加,挿入

push!(a, 13)         #要素と同じ型追加
append!(a, [1,2,3])  #末尾に[1,2,3]を追加
prepend!(a, [1,2,3]) #先頭に[1,2,3]を追加
insert!(a, 3, 15)    #3番目に15を挿入する。(要素数は増える)
vcat(a, [11,12,13])  #aに末尾に[11,12,13]を追加(aは非破壊になる)
hcat(b, [10,11,12])  #bに列を追加
vcat(b, [10,11,12]') #bに行を追加('は転置)
要素数があらかじめわかっているなら、始めに配列を確保して代入したほうが速いっぽい。

置換

a[1] = 15
a[1:4] = [4,4,4,4]
#a[2:5]を削除して[0,0,0]を挿入。戻り値は削除した元のaの[2:5]の値。
splice!(a, 2:5, [0,0,0]) 

削除

shift!(a)  #最初の要素を削除
pop!(a)    #最後の要素を削除

検索

3 in a
findfirst(a, 4) #はじめに4が出る要素番号、なければ0を返す
findin(a, 4)    #4の存在する要素数をArrayに入れて返す)

Lytroが日本に上陸してアクセスが増えたのでLytroのすごいと思うところをのべる。

Lytroが日本にやってくるらしい。
http://dc.watch.impress.co.jp/docs/news/20141105_674603.html

なんか記事のアクセスがすごい。
これは良い機会とばかりにLytroがすごいと思うところを挙げる。1つは文化的内容。もう一つは技術的内容。

1. アナログとデジタルの融合が進んでいる。

上の記事にレンズの断面図がのってある。35mm判換算で30-250でFナンバーは2.0らしい。最近のカメラと比べると広角側が少し足りないと思われるかもしれない。注目すべきはレンズのスペックが高いにも関わらず、枚数12枚で非球面レンズを使っていないことだ。
http://dc.watch.impress.co.jp/img/dcw/docs/674/603/html/17.jpg.html

 確かIllumの発表の時に言及しとったと思うけど、レンズの収差もデジタルで補正できるような最適設定にしているらしい。収差のデジタル補正と言えば富士フイルムとかキヤノンにそういう機能(デコンボリューション)がついていたと思うけど、こういう技術はどちらかというと、まだレンズの最適化で取れない部分をデジタルで補正しているイメージ。しかし、Lytroの場合発想が逆で、デジタルで補正できるところまでレンズの収差を出しているようだ。Ngさんの博士論文にもそういう研究はあったし。もちろん日本の企業でも「倍率色収差」とか「歪曲収差」とかは、デジタルで補正できるようにセンサーとレンズの最適化を行っとるやろうけど、Lytroの最適化は1つ先の次元に行っとる気がする。


2. 技術をもって本気で新しい表現を模索している

Lytroは新しい表現を探している。ホームページのムービーでも言っとるけど、Light Fieldカメラで一瞬を切り取るとどういう表現ができるか模索している。最初に出した製品はセンサーも小さく使い勝手も悪かったが、今回のIllumはクリエイターをターゲットにしている。Light Fieldカメラが生き残るかどうかは、新しい表現方法が見つかるかどうかにかかっているといってもいい気がする。ボケ味を追求する日本のメーカーとは次元が違う。(もちろんどっちがいいとかいう話は抜きにして)
残念ながらぼくにはどういう表現がよいのか全くわかりません。

やっぱりSilicon Valleyの企業は日本とは考え方が違う。
ただLytroのIllumは使い勝手が悪いらしい。やっぱり手に馴染む機器の作り込みには時間がかかるようだ。


2014年11月3日月曜日

Juliaで固定配列

JuliaでArrayを使えばベクトルや行列計算できるけど、2とか3次元の点を表すのに可変長の配列を使うのはどうも無駄が多い気がする。自前で
immutable Vec
    x::Float64
    y::Float64
    z::Float64
end
とかやったけど、演算子とかを定義するのが面倒。

そんな問題を解決してくれるのがImuutableArrays
Julia0.3にはすでに採用されているので、
Using ImuutableArrays
v3 = Vector3{Float64}(1, 2, 3)
v2 = Vector2{Float64}(4, 5)

と書ける。デフォルトでは読み込んだ時に4次元まで作ってくれて、演算子の定義とかもやってくれている。(これがメタプログラミングの力か…)
しかも上で自分で定義したのより速い…

行列は
typealias Mat3d Matrix3x3{Float64}
とすれば使える。






2014年11月2日日曜日

IPython notebook でプレゼンテーション用スライドを作る

プレゼンテーションを作るときに困ることがある。

  1. 手軽に書きたい
  2. Texで数式
  3. 図を入れたい。
  4. プログラムを入れたい。
  5. かっこ良くしたい。
  6. PDFにしたい。
  7. 個人で使う分はフリーで済ませたい。
ざっと思いついたままに書くとこうなる。
いろいろ探していたら、最近IPython notebookからスライドが作成できることが判明した!
しかもreveal.jsを使っているのでいかにもギークっぽいかっこいいスライドが作れる。
こんな感じのやつ→http://lab.hakim.se/reveal-js/#/

まえからreveal.jsは使ってみたかったけど、jsとかHTMLとかよく分からんし、Web上で使えるやつもあったけど、PDFにするのに有料会員にならんといけんとかで挫折した。


スライドの作り方

1. IPython Notebookをインストールする。

IPython NotebookはAnacondaをつかってインストールするのがおすすめ。一式全部はいるし、アップデートも
conda update anaconda
で完了する。

2. 新しいNotebookを開いて設定する。

New Notebookを作成したら以下の図のところをSlideに設定する。

右のタブはSlide、Sub-Slideとか選べる。これはreveal.jsでいう横にスライドか、縦にスライドかを指定できる。いろいろ試してみて。

3. スライドへ変換する

IPython Notebookの編集が終わったら、NotebookからSlideへの変換をする。スライドの作成はここにあるとおりのコマンドを実行するとできる。--postでHTMLサーバーを指定する必要があるらしい。

ipython nbconvert notebook.ipynb --to slides --post serve

IPythonのバージョンが古いと文字がものすごく小さいので最新版にアップデートするのをおすすめ。
※Chromeで開くと数式がレンダリングされなかったのでmathjaxのプラグインをいれた。
ちなみにPDFにするときは
ipython nbconvert notebook.ipynb --to latex --post pdf


問題は記録用にPDFとかHTMLにしたらちょっとダサくなってしまうことか。
まぁよい。


2014年9月25日木曜日

Opencv3.0 ALPHA Python版をMacにソースからインストールしたときのメモ

Mac OSXでPythonはAnacondaでインストールしたPython3.4を使用


  1. OpenCV for Linux/MacをダウンロードしてZipを解凍する
  2. ダウンロードフォルダに移動して、作業用フォルダを作っておく。ここではbuildとした。
  3. CMakeをダウンロードして実行。
  4. ダウンロードしたOpenCVのフォルダとbuildフォルダを指定する。
  5. Searchにpythonと入れてConfigを押す
  6. Configを押すとPython3の設定がいろいろあるので上の図や参考サイトを参考にして埋める。赤い表示がなければとりあえずOK。
  7. Generateを押す。
  8. buildフォルダで
    make -j8
    を実行。失敗するとエラーが出る。Python.hがないとか言われたけど、6の設定したPATHを見直すとうまくいった。
  9. sudo make install
  10. pythonを立ち上げてimport cv2
  11. 今回は「 Library not loaded: libpython3.4m.dylib」とか言われたのでシンボリックリンクを貼った。
  12. sudo ln -/Users/name/anaconda/lib/libpython3.4m.dylib /usr/local/lib

これでとりあえずopencvが動いた。

参考サイト
  1. http://luigolas.com/blog/2014/09/15/install-opencv3-with-python-3-mac-osx/
  2. https://stackoverflow.com/questions/20953273/install-opencv-for-python-3-3

2014年9月21日日曜日

Light fieldの原理をLytroがじきじきに解説

久しぶりにLytroのblogをのぞいたら、Lytroが論文を公表していた。
まだ目を全部に目を通してないけど、見たことない頭があったりした。
http://blog.lytro.com/post/97168596625/modeling-the-light-field-camera

最近はGoogle Cameraでもボケを作れるので静止画でのLytroの優位性はなくなってきている気が。動画だと処理が大変そうだしなー。

2014年9月19日金曜日

Microsoft OCR Library がリリースされた。

https://www.nuget.org/packages/Microsoft.Windows.Ocr/

WindowsRuntimeってWinRTのことで、Windows Metroのアプリをつくるためのものか。
いつか使ってみたいけど、言語がどれもわからん。

2014年9月8日月曜日

Juliaで簡単な回帰

Juliaの練習がてら簡単な線形回帰を行った。

以下のことが分かる練習にはちょうど良いシンプルなコードだ。

  • 乱数をどう発生させるか
  • 計画行列をどう作るか
  • Gadflyを使って2つのプロットを重ねるにはどうするか
  • Gadflyを使って軸の設定をどうするか
これらがわかれば自作する当面の計算にはこまらない気がする…
linreg()はバイアス項を勝手に計算してくれるみたいなので、入れなくてよいみたいだ。

# #Regression Example

using Gadfly
using Distributions

N = 15
xdata = (rand(N) - 0.5) * 2 * pi
σ² = 0.5
function noisy(x, r=rand(Normal(0, σ²))
    return sin(x) + r
end

#X results in Nx9 array
t = noisy(xdata);
X = Float64[x^i for x=xdata, i=[1:9]] #bias is not needed for linreg(). it will interpolate automatically

#regression by linreg
coeffs = linreg(X, t)

#plotting
f = x -> coeffs[1] + sum([coeffs[i]*x^(i-1) for i = [2:10]])

x_reg = linspace(-pi, pi, 100)
y_reg = [f(xi) for xi=x_reg]

plot(layer(x=x_reg, y=y_reg, Geom.line), layer(x=xdata, y=t, Geom.point),
    Coord.Cartesian(xmin=-pi, xmax=pi, ymin=-1.5, ymax=1.5))


2014年8月22日金曜日

グラフ読み取り

PyQtでグラフを読み取れるアプリを作った。
グラフの画像をCtrl+Vでペーストし、その座標をクリックすると画素値からグラフの値を教えてくれる。

最初に右クリックで左下と右上の座標を教える必要がる。
なんでこんなんがいるかって、世の中紙でデータを保存する人がまだいるということである。

(OpenCVは使っていない)

from PyQt4 import QtGui, QtCore, QtCore
import sys
import cv2
import numpy as np

def convert_array_to_qimg(cv_img):
    """
    this is based on https://stackoverflow.com/questions/18676888/how-to-configure-color-when-convert-numpy-array-to-qimage
    """
    height, width, bytesPerComponent = cv_img.shape
    bytesPerLine = bytesPerComponent * width;
    #cv2.cvtColor(cv_img, cv2.CV_BGR2RGB)
    return QtGui.QImage(cv_img.data, width, height, bytesPerLine, QtGui.QImage.Format_RGB888)

class GraphCoordinate:
    def __init__(self):
        pass

    def setBottomPixelCoordinate(self, x1, y1):
        """
        x1 and y1 are pixel
        """
        self.x1 = x1
        self.y1 = y1

    def setBottomGraphCoordinate(self, X1, Y1):
        self.X1 = X1
        self.Y1 = Y1
        print "(X1, Y1) = (%.2f, %.2f)" % (self.X1, self.Y1)


    def setTopPixelCoordinate(self, x2, y2):
        self.x2 = x2
        self.y2 = y2

    def setTopGraphCoordinate(self, X2, Y2):
        self.X2 = X2
        self.Y2 = Y2
        print "(X2, Y2) = (%.2f, %.2f)" % (self.X2, self.Y2)

    def getGraphCoordinate(self, x, y):
        w = x - self.x1 
        h = self.y1 - y
        W = self.x2 - self.x1 
        H = self.y1 - self.y2 
        X = self.X1 + w * (self.X2 - self.X1) / float(W)
        Y = self.Y1 + h * (self.Y2 - self.Y1) / float(H)
        return X, Y 



class MainWindow(QtGui.QMainWindow):
    def __init__(self):
        super(MainWindow, self).__init__()
        self.image_frame = QtGui.QLabel()
        fname = r"sample08.png"
        img = cv2.imread(fname)
        #img = QtGui.QImage(fname)
        qimg = convert_array_to_qimg(img)
        self.image_frame.setPixmap(QtGui.QPixmap.fromImage(qimg))
        self.clip = QtGui.QApplication.clipboard()
        
        self.setCentralWidget(self.image_frame)
        self.move(30,30)

        self.counter = 0
        self.coord = GraphCoordinate()
 
    def keyPressEvent(self, e):

        if (e.modifiers() & QtCore.Qt.ControlModifier):
            #selected = self.table.selectedRanges()
            if e.key() == QtCore.Qt.Key_V:#past
                qimg = self.clip.image()
                self.displayImage(qimg)

            elif e.key() == QtCore.Qt.Key_C: #copy
                pass

    def mousePressEvent(self, e):
        #print e.x(), e.y()
        if e.button() == QtCore.Qt.RightButton:
            if self.counter % 2 == 0:
                text, ok = QtGui.QInputDialog.getText(self, 'Input Dialog', 
                                                        'Enter X1, Y1 (comma separated):')
                X1, Y1 = [float(v) for v in text.split(',')]
                if ok:
                    self.counter += 1
                    self.coord.setBottomPixelCoordinate(e.x(), e.y())
                    self.coord.setBottomGraphCoordinate(X1, Y1)


            else:
                text, ok = QtGui.QInputDialog.getText(self, 'Input Dialog', 
                                                        'Enter X2, Y2 (comma separated):')
                X2, Y2 = [float(v) for v in text.split(',')]
                if ok:
                    self.coord.setTopPixelCoordinate(e.x(), e.y())
                    self.coord.setTopGraphCoordinate(X2, Y2)
                    self.counter = 0

        elif e.button() == QtCore.Qt.LeftButton:
            try:
                x, y = self.coord.getGraphCoordinate(e.x(), e.y())
                print "%.5f %.5f" % (x, y)

            except AttributeError as e:
                print "please set the graph coordinate by right clicks."


    def displayImage(self, qimg):
        self.image_frame.setPixmap(QtGui.QPixmap.fromImage(qimg))

 
def main(args):
    app = QtGui.QApplication(sys.argv)
    form = MainWindow()
    form.show()
     
    sys.exit(app.exec_())
 
if __name__ == "__main__":
    main(sys.argv)


2014年8月10日日曜日

今さら聞けないRichardson-Lucy入門

デコンボリューションがとても魅力的だ。これを使えばレンズの収差やボケ、手ブレによる画像劣化が除去できるという夢の技術だ。

いろんな論文を読むとかならずといっていいほど以下の2つの手法が挙げられている。
1.Wiener Filter
Wikipediaに導出方法が載っている。

2.Richardson-Lucy Deconvolution
こちらもWikipediaに載っている。Wikipediaを読んでもわからなかったので、勉強した内容を少し解説してみる。

わかりやすいように1次元で話を進める。

をピクセルiの真の値(全くボケていないときの値)

をピクセルiで観測された値

とする。全くボケていないときはこの2つが一致する。でも今回はPSF(Point Spread Function)によってボケる。ボケがあるとピクセルjへ行く光は広がりをもってしまうので、ピクセルiへ少しおすそ分けしてしまう。
本来jへいくはずの光のボケでのih量を


とする。Richardson-Lucyではこのp(i,j)は既知だとするイメージは下の図。


ここでZ(i,j)なるものを導入する。Z(i,j)は期待値がP(i,j)λjのポアッソン分布に従いう確率変数。

これはピクセルiに届く光が確率的に変化し、期待値はp(i, j)λjということを意味する。
Wikipekiaにこの量は出てこないけど、元論文などを読むと出てくる。
もともとデコンボリューションはCTとかで研究されていて、CTは原子の崩壊(ガンマ線?)を捉えるので、観測される値が確率的になる、という経緯がある。
(一様な確率で起こる事象がある期間で観測される回数はポアッソン分布になる。詳しくはWikipediaをみてください。)

さて、点が1点だけなら復元も簡単かもしれんけど、周りの点の広がりが重なって観測されるので問題が難しくなる。例えば下の図のように重なるのでピクセルiで値は, 周辺のピクセルからの寄与を足しあわせたのものになる。




ここでzは確率的な値なので、観測結果と期待値から予測することで話を前に進めよう。
ピクセルiで観測された値をつかうと、その内ピクセルjからきた値は、割合で考えると下のような式になる。右辺の分母がiへくる光の総和で、分子がその内jからきたものを表す。





zはもともとjにあったものが周辺のピクセルiに散らばったものとかんがえることができるので、集めてやると真の値λjが分かる。真のzはわからないので推定値を使うとλを推定できてる。


この式にさっきのzの推定値を代入すると

このλjは真の値なので、知り得ない。そこで前回推定した値を使って逐次的に求めることにすれば


というふうになる。

実際に実装するときはWikipediaのようにコンボリューションの関数を使って書いたほうが楽だった。


(2014/11/15 式間違ってましたので修正しました。前に書いてたやつではPSFが対称じゃないとうまくいかなかったorz)

参考
https://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution
http://people.csail.mit.edu/dgreensp/may/deconvlucy.html



Juliaで自作VectorとArray{Float64,1}の速度比較

Juliaの勉強を始めた。
Juliaは数値計算の関数を大量に有しているけど、可変長のArrayをもとにしていて、3x3の行列や3x1ベクトルの単純なものだと自分で定義したほうが高速なのだろうか?という疑問が湧いたので実験してみた。
もしかしたらJuliaは裏でものすごいことをやっていてるかもしれんし。

3x1のベクトルVecを以下のように定義した。
immutable Vec
    x::Float64
    y::Float64
    z::Float64
end
dot(a::Vec, b::Vec) = (a.x*b.x + a.y*b.y + a.z*b.z)
cross(a::Vec, b::Vec) = Vec(a.y*b.z - a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x)

ただ外積やらなんやらを計算する関数を作った。
function cross_loop(N::Int, a::Vec, b::Vec)
    tmp = 0.0
    for i in 1:N
        v = cross(a,b)
        tmp += dot(v,v)
    end 
    tmp/N
end

function cross_loop(N::Int, a::Array{Float64,1}, b::Array{Float64,1})
    tmp = 0.
    for i in 1:N
        v = Base.cross(a,b)
        tmp += Base.dot(v,v) 
    end
    tmp / N
end
Arrayの方はデフォルトで入っているBase.dotやcrossをつかった。

N=10^7くらいにして計算
速度は@timeで計測した。1回目はコンパイルに時間を取られるので2回目の時間を使用。
N = 10^7
a = Vec(1,2,3)
b = Vec(4,5,6)
v = [1.,2.,3.]
u = [4.,5.,6.]

@time cross_loop(N, a, b)
@time cross_loop(N, u, v)

結果はダントツで自作のほうが速い。
elapsed time: 0.013432361 seconds (96 bytes allocated)
elapsed time: 2.079223219 seconds (1280000096 bytes allocated, 43.49% gc time)


ちなみにVecをimmutableでなくtypeで定義すると、それだけでかなり遅くなる。GCが働いているようだ。
elapsed time: 0.419364716 seconds (320000096 bytes allocated, 53.07% gc time)


結論:3次元に特化したベクトル計算では自分で定義してやったほうが速い。




2014年6月20日金曜日

これなら分かる、最尤推定と事後確率最大化(MAP推定)の違い

最尤推定と事後確率最大化。どっちもガウス分布を仮定して対数を取ると2次式になるのでいままでなんとなく「同じようなもん」だと思っていた。

ここでは、観測された値をy、潜在変数をxとして話を進める。

まずはじめに訴えたいのは、どちらの手法も潜在変数xを求めたいということ。


尤度とは、あるモデル(仮説)を仮定したときの観測結果の確率P(y|x)である。xはいろいろな値が考えられるけど、観測結果を最もよく表すようにxを決定すれば現実に即していると言えよう。つまり、P(y|x)をxについて最大化するのが最尤推定。



次に事後確率最大化を考える。
まず単に事後確率といったらP(x|y)なのかP(y|x)なのかよくわからない。どちらも高校でならう事後確率の形式やし。
ここでは、観測がされた後の確率を事後確率と呼ぶ。観測が条件となる確率。つまりP(x|y)のこと。P(x|y)を最大化するの事後確率(MAP推定)だ。では最尤推定とどうちがうのか??

ベイズの定理を使うと
P(x|y) ∝P(y|x)P(x)
右辺は(尤度)*(事前確率)となっている。
ここで(事前確率)を何かの定数だと考えてP(x|y)を最大化を考えると、上の尤度を最大化する場合と同じになっている!
つまり、事後確率最大化は最尤推定に事前確率というバイアスがかかったものだということだ。事後確率は主観とか言われてこれがベイズと呼ばれるアプローチだ。
事後確率の負の対数を取って最大化すると正則化項として出てくる(よく教科書にでてくるパターン)。つまり主観がある種のブレーキとして働くと考えられる。






2014年6月17日火曜日

PyQt: Paste the clipboard image.

The last time I showed how to copy&paste the clipboard text with Excel.
Today I will show you how to deal with a copied clipboard image.

All you have to do is just like below. Not much different from the text example.
 
self.clip = QtGui.QApplication.clipboard()
qimg = self.clip.image()



You might wannat convert QImage to Numpy array, but I couldn't find the very nice solution.
I will just save the image and cv2.imread. This might be slow, but I can accept for now.
qimg.save("tmp.png")
img = cv2.imread("tmp.png")


To covert from numpy array to QImage, I referred this article.
https://stackoverflow.com/questions/18676888/how-to-configure-color-when-convert-numpy-array-to-qimage
def convert_array_to_qimg(cv_img):
    height, width, bytesPerComponent = cv_img.shape
    bytesPerLine = bytesPerComponent * width;
    #cv2.cvtColor(cv_img, cv2.CV_BGR2RGB)
    return QtGui.QImage(cv_img.data, width, height, bytesPerLine, QtGui.QImage.Format_RGB888)


The overall program looks like below.
from PyQt4 import QtGui, QtCore
import sys
import cv2
import numpy as np

def convert_array_to_qimg(cv_img):
    """
    this is based on https://stackoverflow.com/questions/18676888/how-to-configure-color-when-convert-numpy-array-to-qimage
    """
    height, width, bytesPerComponent = cv_img.shape
    bytesPerLine = bytesPerComponent * width;
    #cv2.cvtColor(cv_img, cv2.CV_BGR2RGB)
    return QtGui.QImage(cv_img.data, width, height, bytesPerLine, QtGui.QImage.Format_RGB888)



class MainWindow(QtGui.QMainWindow):
    def __init__(self):
        super(MainWindow, self).__init__()
        self.image_frame = QtGui.QLabel()
        fname = r"sample08.png"
        img = cv2.imread(fname)
        #img = QtGui.QImage(fname)
        qimg = convert_array_to_qimg(img)
        self.image_frame.setPixmap(QtGui.QPixmap.fromImage(qimg))
        self.clip = QtGui.QApplication.clipboard()
        
        self.setCentralWidget(self.image_frame)
        self.move(30,30)
 
    def keyPressEvent(self, e):
        if (e.modifiers() & QtCore.Qt.ControlModifier):
            #selected = self.table.selectedRanges()
                 
            if e.key() == QtCore.Qt.Key_V:#past
                qimg = self.clip.image()
                self.displayImage(qimg)

            elif e.key() == QtCore.Qt.Key_C: #copy
                pass

    def displayImage(self, qimg):
        self.image_frame.setPixmap(QtGui.QPixmap.fromImage(qimg))


 
def main(args):
    app = QtGui.QApplication(sys.argv)
    form = MainWindow()
    form.show()
     
    sys.exit(app.exec_())
 
if __name__ == "__main__":
    main(sys.argv)






2014年6月3日火曜日

scikit-learnを使って数字認識(2) k-NNを使った学習

scikit-learnを使えば簡単な機械学習ならすぐにできる。
今回は特徴ベクトルとして画素値をもちい、それをk-NearestNeighborsにぶちこむ。
k-meansもためしたけど、こちらはあまりうまく行かなかった。
Pythonだと学習した識別器(オブジェクト)をPickleで保存して簡単に取り出せる。

10種類くらいのフォントを学習させて11種類目でためしたが、正答率は100%だった。今後もっといろんな状況下(画素が少ないサンプルとか)で試していってどこまでいけるか確かめたい。

以下メモ

学習部分
import numpy as np
import cv2
from sklearn.neighbors import KNeighborsClassifier
from sampling import convert_to_binary
"""
1. Read sample image and convert to 1d feature vector
2. pass the feature vectors to kmeans clustering maching
3. pickle the result.
"""

def convert_to_feature_vector(src):
    """
    convert source image to a 1d feature vector
    """
    N =10
    im = cv2.GaussianBlur(src,(5,5),0)
    im = cv2.resize(src, (N, N))
    im = im.reshape(N**2)
    #im = np.array(im>124, dtype=np.int8) #convert to 0 and 1
    return im


if __name__ == "__main__":
    #read jpgs, and resize them as a vector

    dirname = [str(i) for i in range(10)]
    dirname += ["dot"]
    dirname += ["bar"]

    X = [] #sample data
    Y = [] #label
    for dn in dirname:

        #input label(0-11)
        if dn == "bar":
            label = 11
        elif dn == "dot":
            label = 10
        else:
            label = int(dn)


        fnames = os.listdir(dn)
        fnames.remove(".DS_Store")

        for fn in fnames:
            #print os.path.join(dn, fn)
            im = cv2.imread(os.path.join(dn, fn))
            im = cv2.cvtColor(im, cv2.COLOR_BGR2GRAY)
            im = convert_to_feature_vector(im)
            #print im.reshape(10,10)
            X.append(im)
            Y.append(label)

        #convert to np.array
    X = np.array(X)

    #
    knn = KNeighborsClassifier(n_neighbors=5)
    knn.fit(X, Y)

    #save the classifier as pickle
    import pickle
    with open("knn_trained.dump", "w") as f:
        pickle.dump(knn, f)



予測部分
#coding: utf-8

from sklearn import cluster
import pickle
import cv2
from training_knn import convert_to_feature_vector


if __name__ == "__main__":
    with open("knn_trained.dump") as f:
        knn = pickle.load(f)
        print knn

    import os
    dirname = "test_data"
    fnames = os.listdir(dirname)
    try:
        fnames.remove(".DS_Store")
    except ValueError as e:
        print e


    for fn in fnames:
        im = cv2.imread(os.path.join(dirname, fn))
        im = cv2.cvtColor(im, cv2.COLOR_BGR2GRAY)
        x = convert_to_feature_vector(im)
        y = knn.predict([x])
        print fn, y



scikit-learnを使って数字認識(1) サンプル取得

OCRはフリーのものもあるけど、数値データを起こしたいときには無駄に高機能になってしう。それだけでなくご認識も増える。

手書きでない数値(0-9, "-", ".")を認識するだけならそんなに苦労しないのでは?と思ったのでちょっと実験してみた。

フツウの人は手書き文字認識に興味があるやろけど、
調べてみるとStackoverflowに以下のような投稿があって結構盛り上がっている。
https://stackoverflow.com/questions/9413216/simple-digit-recognition-ocr-in-opencv-python

今回は

  1. 上のリンクの方法で数値を切り取り
  2. scikit-learnを使って学習
  3. 学習データを元に新しい予測
ということを行う。

まず、以下のように切り取られた画像から数値を取得したい。


イメージしている成果物は、WebやPDF上をマウスで矩形選択するとそこの数値をテキストデータに変換するツールだ。
ほぼ上のサイトと同じことをしているが、一つだけ違うところfindContoursで輪郭の取得オプションだ。opencvでは輪郭の取得のみならず、それら輪郭の関係まで返してくれる。今回はcv2.RETR_CCOMPを指定した。(一番外側の輪郭とするとなぜか図の枠が検出されたので)輪郭についてはこのページがわかりやすかった。輪郭の階層を指定することで輪郭の大きさによる分類など無駄な作業が省ける。
下のプログラムを実行するとこうなる。


if __name__ ==  "__main__":
    sample_dir = "fonts"
    sample_file = "sample03.png"
    im = cv2.imread(os.path.join(sample_dir, sample_file))
    im_copy = im.copy()
    gray = cv2.cvtColor(im,cv2.COLOR_BGR2GRAY)
    blur = cv2.GaussianBlur(gray,(5,5),0)
    thresh = cv2.adaptiveThreshold( src=blur, 
                                    maxValue=255, 
                                    adaptiveMethod=cv2.ADAPTIVE_THRESH_MEAN_C,#cv2.ADAPTIVE_THRESH_GAUSSIAN_C,
                                    thresholdType=cv2.THRESH_BINARY,
                                    blockSize=5,
                                    C=3)
    thresh_copy = thresh.copy()  #thresh are destroyed when findCountours

    cv2.imshow('threshold', thresh)
    cv2.waitKey(0)


    contours,hierarchy = cv2.findContours(thresh_copy, cv2.RETR_CCOMP, cv2.CHAIN_APPROX_SIMPLE)
    print "%d contours are found" % len(contours)


    i = 0
    for con, hie in zip(contours, hierarchy[0]):
        if hie[3] != -1: #2:first child, 3:parent
            (x,y,w,h) = cv2.boundingRect(con)
            #roi = thresh[y:y+h,x:x+w]
            roi = im[y:y+h, x:x+w]

            cv2.rectangle(im_copy,(x,y),(x+w,y+h),(0,0,255),2)
            cv2.imshow('rect',im_copy)
            #cv2.imshow("threshold", roi)
            
            key = cv2.waitKey(0)
            if key == 27: # (escape to quit)
                sys.exit()

            c = chr(key)

            if c == ".":
                c = "dot"
            elif c == "-":
                c = "bar"
            elif c == "n": #negative
                c = "negative"

            file_head = sample_file.replace(".png", "")
            if c == "negative":
                fname = os.path.join(c, file_head + "_" + c + ("%02d.png"%i))
                i+=1
            else:
                fname = os.path.join(c, file_head + "_" + c + ".png")

            #write ROI 
            print fname
            #cv2.imwrite(fname, roi)









2014年5月28日水曜日

とりあえずBoost.Pythonを使う(5) return_internal_referenceで==をしてみると…

Boost.Pythonで保有しているオブジェクトをgetterなどで返すときはreturn_internal_referenceを指定すると便利だ。
class_<Surface>("Surface", init<>())
        .def("__str__", &Surface::toString)
        .add_property("shape", 
                      make_function(&surface_get_shape, return_internal_reference<>()),
                      make_function(&surface_set_shape, with_custodian_and_ward_postcall<1, 2>())
                      )

こうするとオブジェクトを取り出して値を変更すると内部の値も変わっている。
唐突だが、以下の例では途中で片方の値だけ変更したけど、もう片方も変わっている。
shape = eva.Sphere(13,2)
self.surf.shape = shape

a = self.surf.shape
b = self.surf.shape
print(a,b)

a.r = 1000
print(a,b)

出力結果:
Sphere r: 13, erad: 2 Sphere r: 13, erad: 2
Sphere r: 1000, erad: 2 Sphere r: 1000, erad: 2

これは想定内の働き。でもここで
a == b
とすると、Falseになってしまう。なぜだ。
中身をみるとどうやらアドレスが別らしい。BoostPython側で何か処理をしてくれていてこうなったのだろう。


こうなるとオブジェクトが等しいという判断は、オブジェクトのアドレスでなく、中身の値で判断させるようにしたほうが良さそうだ。

そこで必ず==を実装するようにすべく、
C++のVirutualクラスに==や!=をオーバーライドさせるように書く。
引数がネックになってどう書いていいかわからなかった。
書き方を調べると以下のような記事が見つかった。
https://groups.google.com/forum/#!topic/comp.lang.c++.moderated/VAm49RqRico


メモがてら抜粋させていただくと
bool Derived::operator==(const Base &rhs) const
{
    const Derived *pRhs = dynamic_cast<const Derived *>(&rhs);
    if (pRhs == NULL) { // rhs is of another type derived from Base
         return false;
     } else {
          // Do work to compare two Derived objects
     }
}
違うクラスが来るとfalseになるような仕様になっている。
コンパイルも通ったし、めでたしめでたし。

ただ、質問の解答は下に続いていて、
そもそもPure Virtualクラスに==をオーバライドさせるとか設計が悪すぎるとか書いてある。んーどうすれば一番良いのだろうか…





2014年5月25日日曜日

とりあえずBoost.Pythonを使う(4) pythonのリストと関数のオーバーライド

前にも書いた気がするけど、Boost PythonでC++のベクトルを返すときにvector_index_suiteを使えば受け渡しが簡単だ。だけど、Pythonのリストとして渡したり、返したりしたほうがPython側からすればシームレスに扱えてよい。そんな時はboost::python::listを使う。
例えば
void Hoge::setX(std::vector<double> &x) {
    _x = x;
}

というC++のメソッドがあった場合Pythonでは
>>> hoge.setX([1,2,3])
としたい。

そんなときはboost::python::listを引数にして、boost::python::extrac<>で値を抽出すればよい。
void Hoge::setX(boost::python::list &x) {
    std::vector<double> tmp;
    int len = boost::python::len(x);
    for(int i=0; i<len; i++) {
        tmp.push_back(boost::python::extract<double>(x[i]));
    }
    _x = tmp;
}

こうやってしまうと元の関数を書きなおしてしまうでのC++を使う場合には使いにくい。
同じ名前のメソッドを作ってその内1つをPythonのWrapperに登録する。
static void hoge_set_x(Hoge &h, boost::python::list &x){hoge.setX(x);}

boost::python::class_<Hoge>("Hoge")
               .def("setX", &hoge_set_x)
               ;

ただ、このやり方だと、Wrapperを2つも各のでちょっと周りくどい。しかもlistとvectorをいちいちコピーしないといけない。index_suiteを使った良いやり方がありそうな気がするけど、std::vector<std::vector<> >でエラーがでたので、今回は諦めて上のようにしてみた。

containerのラッパーを書いている以下のようなサイトも見つかった。
https://wiki.python.org/moin/boost.python/StlContainers
大変なので今後の参考としておこう。





2014年5月18日日曜日

とりあえずBoost.Pythonを使う(3) Abstract Class

Boost.Pythonを使ってC++の継承関係をラップすることもできる。
親クラスがvirtual=0なメソッドを持っているときは、(Abstract ClassとかInterfaceとか)は以下のようにすればよい。
class BaseShape {
public:
    virtual ~BaseShape(){}
    virtual std::string toString() const = 0;
};

class Sphere : public BaseShape {
public:
    Sphere(double r);
    ~Sphere();
    std::string toString() const;

private:
    double _r;
};

ラッパー側は以下のようにする。
    class_<BaseShape, boost::noncopyable>("BaseShape", no_init)
        .def("__str__", &BaseShape::toString)
        ;
    class_<Sphere, bases<BaseShape> > ("Sphere", init<double>())
        .def("__str__", &Sphere::toString)
        ;

no_initを記述してなくて、コンパイル時ずっとabstract classがどうとかいうエラーが出ていた。

Boost.PythonをPython3で使う。

HomebrewでインストールされるBoost.pythonは「python」コマンドで呼び出されるpythonを使ってコンパイルさえるのでデフォルトでは2になっている。

でも今後はPython3で開発したいので、Boost.PythonもPython3対応にしてみる。
HomebrewではPython3の指定の仕方がわからないので自分でビルドする。ビルドとなると一歩引いてしまう自分がいたが、簡単だった。

まずboostをダウンロードする。
そのフォルダに移動して
./bootstrap.sh --with-libraries=python --with-python=python3 --with-python-version=3.3

とすれば設定ファイルproject-config.jamが出力される。--with-librariesで欲しいものだけコンパイルできるようだ。出力されたファイルをみてPython3になっていればよい。

それで
./b2

とすればコンパイルが始まる。エラーが結構でたきがするけど、気にしないでおこう。
最後にメッセージが出る。

The following directory should be added to compiler include paths:

    /Users/yourname/Downloads/boost_1_55_0

The following directory should be added to linker library paths:

    /Users/yourname/Downloads/boost_1_55_0/stage/lib


指定されたフォルダを見て、libboost_python3.aがあればOK.
コンパイル時にそれをリンクすればPython3とC++の連携ができる。

2014年5月14日水曜日

とりあえずBoost.Pythonを使う(2) dictを使う

Boost.Pythonを使ってC++とPython間でdictionaryのやり取りを行う。
まずはじめに出会う例はmap_index_suiteを使う例だが、これはpythonで使う場合に.keys()などのメソッドが使えないので困る。
なのでboost::python::dictを直接やり取りするようなメソッドを準備する。

まずPythonのdictionaryをC++のmapにセットする場合のクラス・メソッドは例えば以下のようになる。

#include <boost/python.hpp
void Hoge::set_dict(const boost::python::dict &pydict)
{
    std::map<std::string, double> newmap;

    int len = boost::python::len(pydict);
    boost::python::list keys = pydict.keys();

    for(int i=0; i < len; i++) {
        std::string k = boost::python::extract<std::string>(keys[i]);
        newmap[k] = boost::python::extract<double>(pydict[k]);
    }
    _dict = newmap;
}

逆にC++のマップをPythonのdictionaryとして返す場合は以下のようにする。
boost::python::dict Hoge::get_dict() const
{
    boost::python::dict pydict;

    std::map<std::string, double>::const_iterator it;
    for(it = _dict.begin(); it != _dict.end(); it++)   
          pydict[it->first]=it->second;        
    
    return pydict; 
}




2014年5月13日火曜日

とりあえずBoost.Pythonを使う(1)

Boost.PythonをつかってC++のクラスをPythonへ渡す方法。
この例では、

  1. メンバ変数
  2. メソッド
  3. コンストラクタ
  4. オペレータ
  5. Pythonリストの受け渡し

が分かる仕組みになっている。
例えば3次元ベクトルクラスVector3dをC++で書いたとする。Pythonで扱うとしたら以下のようにしたい。

 v1 = vector3d(1,2,3)   #constructor
 v2 = vector3d([1,2,3]) #constructor with list
 v1.x                   #member
 v1.norm()              #method
 v1+v2                  #operator


ラッパーは以下のように作ってみた。
Pythonのlistとやり取りするクラスが始めから備わっているようだ。
extract<>についてはhttps://wiki.python.org/moin/boost.python/extract


#include <math.h>
#include <boost/python.hpp>

class Vector3d {
public:
    double x;
    double y;
    double z;

    //constructor
    inline Vector3d(double x, double y, double z):
            x(x), y(y), z(z) {}

    //constructor with python linst
    inline Vector3d(boost::python::list pl)
    {
        x = boost::python::extract<double>(pl[0]);
        y = boost::python::extract<double>(pl[1]);
        z = boost::python::extract<double>(pl[2]);

    }
    
    inline Vector3d operator+(const Vector3d &v) const
    {
        return Vector3d(x+v.x, y+v.y, z+v.z);
    }

    inline Vector3d& operator+=(const Vector3d &v)
    {
        x += v.x;
        y += v.y;
        z += v.z;
        return *this;
    }

    //return l2 norm
    inline double norm() const
    {
        return sqrt(x*x + y*y + z*z);
    }

    //to python list
    inline boost::python::list toList() const
    {
        boost::python::list pl;
        pl.append(x);
        pl.append(y);
        pl.append(z);
        return pl;
    }

};



BOOST_PYTHON_MODULE(hoge)
{
 class_<Vector3d>("vector3d", init<double, double, double>())
  .def(init<boost::python::list>())
  .def_readwrite("x", &Vector3d::x) //member variables
  .def_readwrite("y", &Vector3d::y)
  .def_readwrite("z", &Vector3d::z)
  .def(self + self)    //operator
  .def(self += self)
  .def("norm",     &Vector3d::norm) //method
  .def("tolist", &Vector3d::toList) //to python list
  ;
}






2014年5月11日日曜日

とりあえずBoost.Pythonを使う(0)  Makefile for Homebrewer

Boost PythonをOSXで使ってC++モジュールをPythonから使うことができた。

http://d.hatena.ne.jp/moriyoshi/20091214/1260779899を参考に始めた。

上のサイトのコンパイルオプションだと、コンパイルはとおるけど、エラーが出てしまう。
>>> import basic
>>> a = basic.accumulator()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: __init__() should return None, not 'NoneType'

BoostPythonをビルドしたときに使用したPythonを使用しないとコンパイルは通るけど、使用時にエラーが起きるようだ。なので以下のようなMakefileを作ったらうまくいった。

#makefile for boost python
#PYTHON_ROOT = /usr/local/Cellar/python3/3.3.5/Frameworks/Python.framework/Versions/3.3
PYTHON_ROOT = /usr/local/Cellar/python/2.7.6/Frameworks/Python.framework/Versions/2.7

# location of the Boost Python include files and library
 
BOOST_INC = /usr/local/Cellar/boost/1.55.0_1/include
BOOST_LIB = /usr/local/Cellar/boost/1.55.0_1/lib
 
# compile mesh classes
TARGET = basic
 
$(TARGET).so: $(TARGET).o
    g++ -bundle -Wl $(TARGET).o -lboost_python -L$(BOOST_LIB) -L$(PYTHON_ROOT)/lib -lpython2.7 -o $(TARGET).so
 
$(TARGET).o: $(TARGET).cpp
    g++ -I$(PYTHON_ROOT)/include/python2.7 -I$(BOOST_INC) -fPIC -c $(TARGET).cpp



brew でboost pythonをインストールする。

PythonからC++を呼び出せるBoost.Python。Ubuntuではすんなり出来たけど、Macではなかなかうまく行かない。コンパイルは通ってもimportすると、Fatal Python errorとか言われる。
結局
brew uninstall boost
brew install boost --with-python
でimportできるようになった。--with-pythonがないとコンパイル時の-lboost_pythonフラグに対して「boost_pythonなんてライブラリないよー」と言われる。

python3でやってみたいけど、今回はつかれたのでしばらくは2でいこう。

2014年5月4日日曜日

API Design for C++ を読んで学んだことをまとめておく。

順不同。詳しいことは本を買って読んで。
  • Virtual ClassのデストラクタはVirtualで宣言しておく。
    • 子クラスでのみ宣言すると、親クラスの型で変数を定義したときに困る。(P87)
    • メソッドとかは親クラスの型で戻り値を指定しても、子クラスを返すことができる。
  • クラスでメモリを扱う場合はCopy Constructor, Assignment Operator, Destructorの3つをセットで定義すべし。さもないとコンパイラが自動でShallow copyを行う。(P176)
  • getter/setterメソッドは何度も呼び出されるなら、inlineで宣言すればオーバーヘッドがなくなる。でもinlineの中にforとかが会った場合は、勝手にコンパイラがinlineでなくすときがある。小さくて、シンプルなものはinlineするとよい。

boost::tupleを使って2つの戻り値を返すときのメモ

boost::tupleを使って2つの戻り値を返すときのメモ。
作り方はboost::make_tuple(type,1 type2)
取り出し方はinstance.get<i>()

#include <iostream>
#include "eigen-eigen-6b38706d90a9/Eigen/Dense"
#include <boost/tuple/tuple.hpp>

using namespace Eigen;
using namespace std;

boost::tuple<Vector3d, Vector3d> func() {
    Vector3d v1 = Vector3d(1,2,3);
    Vector3d v2 = Vector3d(4,5,6);
    return boost::make_tuple(v1, v2);
}

int main() {
    boost::tuple<Vector3d, Vector3d> vs = func();

    cout << vs.get<0>() << endl;
    cout << vs.get<1>() << endl;
    return 0;
}



2014年5月2日金曜日

Makefileのメモ


以下はMakefileを学習したときのメモ
  • http://shin.hateblo.jp/entry/2012/05/26/231036を参考にした。
  • -oで実行フィアル名の指定。
  • -c オプションでコンパイルのみ実行。リンクはしない。
  • -gでデバッグ情報付加。-Wallでwarningをすべて表示する。
  • ターゲットを分けて書くことで、一部を更新したときに全部をコンパイルしなくて済む。
  • clean: でファイルを消去したりできるけど、このときcleanというファイル存在するとややこしくなるので、.PHONY: clean
  • headerファイルの更新を反映させるために依存ファイルに.hを加える。
  • $@でターゲット名を表す
  • $< で依存ファイルの先頭を表す。
  • $^で依存ファイルのリスト
  • suffix ruleを使えば、拡張子の同じものはまとめて記述できる。
最終的に以下のようになった。といっても上のリンクと同じ。

# Makefile
program = hello
objs = hello.o
CXX = g++
CFLAGS = -g -Wall

#.PHONY: all
#all: $(program)

#primary target
$(program): $(objs) 
    $(CXX) -o $@ $^

#suffix rule
.SUFFIXES: .cpp .o

.cpp .o:
    $(CXX) $(CFLAGS) -c $<

#headers
#hello.o: hello.h

.PHONY: clean
clean:
    rm -f $(program) $(objs)

2014年4月20日日曜日

Raspberry Pi と MPU-6050を使ってデータ取得(C++)

CやC++でMPU-6050のデータを取得したい場合は
http://www.raspberrypi.org/forums/viewtopic.php?t=22266
のプログラムを使えばよい。

RaspberryPiのBを使っている場合はリンク先にも書いてあるようにI2Cdev.cppの中の/dev/i2c-0を/dev/i2c-1に書き換えなければならない。

3D fireframeを動かすでもはみていないけど、クォータニオンを使っているようだ。勉強してないので分からんな。

2014年4月19日土曜日

Raspberry Pi と MPU-6050で角度、加速度データを取得する。

ポータブルな角度データ取得装置を作りたかったので、RaspberryPiとMPU-6050で作ってみた。
MPU-6050はSparkfunとかでは$40とかいう値段やけど、アマゾンだと480円であった。

解説はここに乗っているのでAruduinoを使う場合は参考になりそうだ。
今回の用途はデータを記録したいので、記録容量の観点からRaspberryPiを利用してみる。

まずはこのサイトを見ながら接続と動作を確認する。
Pythonでいけるんですねー。
http://blog.bitify.co.uk/2013/11/interfacing-raspberry-pi-and-mpu-6050.html
http://blog.bitify.co.uk/2013/11/reading-data-from-mpu-6050-on-raspberry.html

今売っているのはRevision2なので注意。
この通りにやるとちゃんと動いていそうだ。


Androcitiには5V動作もOKとか書いてあったのではじめはPiの5Vに接続しとったけど、
データシートをみるとそんなことは書いていない…一応まだ動いている。

Pythonスケールするところは、Rawデータを変換する係数らしい。データシートから値が分かる。
リンクのプログラムでは加速度データを整数で割り算している。これは精度の問題で小数点以下が意味を成さないということなのだろう、多分。止めて測っても下一桁は変わっとるし。




2014年4月16日水曜日

エピポーラ幾何

OpenCVでステレオカメラのキャリブレーションを行うと、引数にEssential MatrixやらFundamental Matrixやら出てきてわからなかった。
英語の文献を読むとわかりやすかった。分かれば大したことないかもしれんけど、ぼくみたいな誰かの参考になるようここにメモっておこう。

参考文献: http://www.robots.ox.ac.uk/~vgg/hzbook/hzbook1/HZepipolar.pdf



ではエピポーラ幾何について話す。
下の図のようにカメラが2つあって、どちらにも点Pが写っているとする。

左のカメラから見た点が、右のカメラのどの座標に写っているかが分かれば、
△OlOrPからPの深さが分かる。

しかし、実際はその対応関係がわからない(下図)。
対応付けをするためには、画素を調べて同じ点を探す必要がある。



しかし、2次元画像のすべての画素の対応を調べる必要はない。
2つのカメラの位置関係がわかっていれば、
その探索は1次元(エピポーラライン上)で済む。

ではその理由を数式を用いて探っていく。
△OlOrPの形は定かではないが、
,


は同じ平面上にあるという制約(Epipolar constraint)がある。
よって以下の式が成り立つ。



これを左のカメラの座標で書き直すと



t、Rは右の座標から左の座標へのtranslation vectorとrotation matrix。
ここで外積はベクトル→ベクトルの変換なので等価な行列表記にすると以下のようになる。





このEがEssential Matrixで、Rとt(カメラの外部パラメータ)から成り立っている。

これはフィルム上の点でも成り立つ。つまり



などとフイルム上の座標を表記したとき、



が成り立つ。


Fundamental Matrixはにカメラの内部パラメータも押し込んだ形のものだ。
フィルム座標からピクセル座標への変換はAffine transform matrix以下のように表せる。



これを上の式に代入すると、





Fはカメラの外部パラメータ、内部パラメータをひっくるめたものであることがわかる。

で、ここから何がわかるかだ。

ピクセル座標での直線を考えるよう。直線の方程式は以下のようになる。




とおくと上のFundamental Matrixの式と同じになる。がライン上にあることを示している。はカメラの構成によってのみ決まるので、どちらかのピクセル座標を与えれば、どのライン上にもう片方のピクセルが存在するか決めることができる。


OpenCVに実装されているので、実際に計算することはないと思うけど、意味は知っておきたいですね。