OpenCV2でSIFT, SURFによる画像の対応付け
大学の講義で画像間の対応付けプログラムを作成するという課題が出されて、無事提出できたのでメモしておこうと思います。このような画像工学のプログラムは初めてで、特にMacOSX + Xcode4 + OpenCV2のサンプルが少なくて色々苦労しました。簡単ではありますがOpenCVのインストールから2画像の対応付け画像を生成するところまでを説明したいと思います。
執筆時のバージョン
OSX 10.7.4 + Xcode4.1 + OpenCV2.4.1
準備
Xcodeはインストール済みであるとして、まずはOpenCVをインストールします。
私はMacPortsを使用しましたが各自で好きな方法でインストールしてください。MacPortsなら
port install opencv
で楽チンにインストールできます
インストールが終わったらXcodeでプロジェクトを作成して、Command Line ToolでC++を選択します。
プロジェクトの作成が終わったらBuild Settingsタブを選択して、
- Valid Architecturesを[i386, x86_64]を[x86_64]に変更
- Header Search Pathsに/opt/local/include/を追加
- Library Search Pathsに/opt/local/lib/を追加
次はOpenCVのライブラリの設定。Build Settingsの隣にあるBuild Phasesタブを選択し、Link Binary With Librariesという項目の＋マークでOpenCVの必要なライブラリ(.dylib)を追加します。私の環境ではopencvとフィルタリングすると一つも見つからなかったので(MacPortsだと/opt/local/libにインストールされるから？)Add Other…から/opt/local/libにあるライブラリを選択しました。今回のプログラムで必要なのは
- libopencv_core.2.4.1.dylib
- libopencv_features2d.2.4.1.dylib
- libopencv_highgui.2.4.1.dylib
- libopencv_imgproc.2.4.1.dylib
- libopencv_nonfree.2.4.1.dylib
（多分）この５つです。バージョン番号は各自でインストールしたバージョンに読み替えてください
プログラム
設定が終わったら後はプログラムを書くだけです。画像間の対応付けには特徴量が必要ですが、これはOpenCVで実装されているSIFTとSURFを使用しました。SIFT、SURFがどのような理論で画像間の対応を取るのかについては中部大学の藤吉教授の論文やスライドが日本語で分かりやすいのでオススメです。リンク
SIFTやSURFで特徴点を求めて、対応付けを行なって対応線を描画するということをマジメにやるとかなり大変なはずですが、ほとんどの機能がOpenCVに実装されているのでそれを使うとかなり少ない行数で実現することができます。
#include <iostream> #include <opencv2/highgui/highgui.hpp> #include <opencv2/nonfree/nonfree.hpp> #include <opencv2/legacy/legacy.hpp> //BruteForceMatcheに必要。opencv2.4で移動した？ int main(int argc, char *argv[]) { //画像読み込み cv::Mat colorImg1 = cv::imread("box1.jpg"); cv::Mat colorImg2 = cv::imread("box3.jpg"); if(colorImg1.empty() || colorImg2.empty()){ std::cout << "No Image" << std::endl; return -1; } //特徴点抽出用のグレー画像用意 cv::Mat grayImg1, grayImg2; cv::cvtColor(colorImg1, grayImg1, CV_BGR2GRAY); cv::normalize(grayImg1, grayImg1, 0, 255, cv::NORM_MINMAX); cv::cvtColor(colorImg2, grayImg2, CV_BGR2GRAY); cv::normalize(grayImg2, grayImg2, 0, 255, cv::NORM_MINMAX); //SIFT // cv::SiftFeatureDetector detector; // cv::SiftDescriptorExtractor extractor; //SURF cv::SurfFeatureDetector detector(1000); cv::SurfDescriptorExtractor extractor; //画像から特徴点を検出 std::vector<cv::KeyPoint> keypoints1; detector.detect(grayImg1, keypoints1); std::vector<cv::KeyPoint> keypoints2; detector.detect(grayImg2, keypoints2); //画像の特徴点における特徴量を抽出 cv::Mat descriptors1; extractor.compute(grayImg1, keypoints1, descriptors1); cv::Mat descriptors2; extractor.compute(grayImg2, keypoints2, descriptors2); //特徴点の対応付け std::vector<cv::DMatch> matches; cv::BruteForceMatcher<cv::L2<float> > matcher; matcher.match(descriptors1, descriptors2, matches); //ソートしたn番目までの対応線を表示させる。nth_elementは要素を基準要素よりも手前に移動させるある種のソート int N=50; nth_element(matches.begin(), matches.begin()+N-1, matches.end()); matches.erase(matches.begin()+N, matches.end()); //対応づけされた画像の用意 cv::Mat matchedImg; cv::drawMatches(colorImg1, keypoints1, colorImg2, keypoints2, matches, matchedImg); /// 画像を表示するウィンドウの名前，プロパティ // CV_WINDOW_AUTOSIZE : ウィンドウサイズを画像サイズに合わせる // CV_WINDOW_FREERATIO : ウィンドウのアスペクト比を固定しない cv::namedWindow("image", CV_WINDOW_AUTOSIZE|CV_WINDOW_FREERATIO); // ウィンドウ名でウィンドウを指定して，そこに画像を描画 cv::imshow("image", matchedImg); // キー入力を（無限に）待つ cv::waitKey(0); return 0; }
OpenCV2.0からC++のインターフェースが追加されたので全てC++で書いてみました。
注意点はOpenCV2.4からSIFTとSURFがfeature2dからnonfreeに移った点です。これに気が付くのに時間がかかってサンプルや他の方のプログラムがコンパイルできないという状況に延々とハマっていました・・・。
その他に手軽に特徴点間の距離を計算して対応を付けを行うBruteForceMatcherのヘッダーがいつからかlegacyに移ったらしいので、ここも注意です。自分は上のnonfreeと合わせてずっと悩んでいましたが、結局stackoverflowで検索して解決に至りました(~_~)
プログラムが何をやっているかは行数も多くないし、コメントを付けてるのでそれほど難しくなく理解できると思います。SIFTとSURFの切り替えはコメントアウトの切り替えで行うことができます。対応線の数はint N=50を変えることで調節できます。類似度が高い特徴点同士から上位N本の対応線を引くので、画像ごとに見やすい数に調節してください。
実行するとこんな感じになります
上：SIFT 下：SURF
SIFTとSURFは回転と拡大縮小に強いので、ある物体が他の画像ではどこに写っているかをこのように見つけることができます。SIFTとSURFの違いは詳しく分かっていないのですが、SURFの方がSIFTを簡略化して計算量を少なくしたもののようです。でも画像で見た感じでは精度にそれほど変わりがないように見えます
上：SIFT 下：SURF
こちらはカメラを平行移動させて、半分ほど共通部分が存在する画像間で対応を取った場合です。両方に写っている部分の対応はかなりいい感じに取れています。いわゆるパノラマ画像はこのようにして画像間の対応を求めてから合成して作られたりするようです.
画像工学のことはほとんど知らなかったので色々と勉強になりました。特に画像を特徴ベクトルで表現する事で今回のプログラムのように画像間の対応付けをしたり、類似画像を探したりすることができることを初めて知りましたが、これって自然言語処理の文書をベクトルで表現するという手法にそっくりですね。音声、自然言語処理、画像、機械学習の分野は手法が似ていると聞いたことがありますが納得。
NumPyとCythonを組み合わせると爆速！
前回の記事の最後にcythonとnumpyを組み合わせても速くならなかったと書いてしまったのですが、@frontier45 さんから公式のチュートリアルをちゃんと読みましょう。と教えていただいたので、自分の勉強がてらブログにも書いておきます。
使用するコードは何のひねりもなく公式のチュートリアルから丸パクリです。
注意：以下のベンチマークはMacOS 10.7.3 MacBook Core2 Duo 2.26GHzでPythonとC/C++のコンパイラはMac標準のPython2.7.1, llvm-gccを使用しています。
Python/Cythonはipythonのtimeitを利用して実行時間を測定しています。
なお、以下の記録はあくまで私の環境、私の実装での記録なので比較の結果は正しいというものではないです。公式のチュートリアルでも実行時間が書かれていますので、そちらの方もご覧下さい
1. コードの変更無し
まずは普通のpythonとただcythonでコンパイルしただけの結果を比較します。前回の結果同様、コードを何も変更しなくても少しは速くなります。
convolve_py.py, convolve_cy.pyx
from __future__ import division import numpy as np def naive_convolve(f, g): if g.shape[0] % 2 != 1 or g.shape[1] % 2 != 1: raise ValueError("Only odd dimensions on filter supported") vmax = f.shape[0] wmax = f.shape[1] smax = g.shape[0] tmax = g.shape[1] smid = smax // 2 tmid = tmax // 2 xmax = vmax + 2*smid ymax = wmax + 2*tmid h = np.zeros([xmax, ymax], dtype=f.dtype) #pythonではxrangeの方が高速なのでここだけ元のコードから変更 #cythonならばどちらでも変わらないはず for x in xrange(xmax): for y in xrange(ymax): s_from = max(smid - x, -smid) s_to = min((xmax - x) - smid, smid + 1) t_from = max(tmid - y, -tmid) t_to = min((ymax - y) - tmid, tmid + 1) value = 0 for s in xrange(s_from, s_to): for t in xrange(t_from, t_to): v = x - smid + s w = y - tmid + t value += g[smid - s, tmid - t] * f[v, w] h[x, y] = value return h
実行結果
In [1]: import numpy as np In [2]: N = 100 In [3]: f = np.arange(N*N, dtype=np.int).reshape((N,N)) In [4]: g = np.arange(81, dtype=np.int).reshape((9,9)) In [5]: import convolve_py In [6]: timeit -n2 -r3 convolve_py.naive_convolve(f,g) 2 loops, best of 3: 2.13 s per loop In [7]: import convolve_cy In [8]: timeit -n2 -r3 convolve_cy.naive_convolve(f,g) 2 loops, best of 3: 1.74 s per loop
2. 型情報の追加
次はcdefによって型情報を加えます。
変更する部分を解説していきます。
#numpyをcythonから扱えるようにするためにcimportを使用します cimport numpy as np #numpyの配列の型情報を扱いやすくするための代入？ DTYPE = np.int #Cにおける#defineに相当 ctypedef np.int_t DTYPE_t
numpy配列の型をCから扱うときは通常の型の名前に"_t"を付けるそうで、numpy.intならnumpy.int_t, numpy.int32ならnumpy.int32_tという具合です。このプログラムでは結果を格納するnumpyの配列 hがnumpy.intで、そこに代入するための計算用変数 valueがnumpy.int_tになっています。
#引数の型がnp.ndarrayであることを明示すると呼び出し時に型チェックが行われたり #f.shape[0]のようなアクセスがより効果的になるらしい def naive_convolve(np.ndarray f, np.ndarray g) #変数はなるべく全てcdefにしておく cdef int vmax = f.shape[0] … cdef int x, y, s, t, v, w ...
とくにforループと配列の要素番号に使用するものをcdefしていないとpythonの遅いループのままなので、忘れずにcdefしておきましょう
実行結果
In [9]: import convolve_ndarray In [10]: timeit -n2 -r3 convolve_ndarray.naive_convolve(f,g) 2 loops, best of 3: 2.68 s per loop
えー！？元のpythonよりも遅いってどういうこと・・・
チュートリアルでの結果と異なってしまい、これについては結局よく分からなかったのですが、ndarray以外の変数だけをcdefにしてやるとちゃんと高速になります
In [11]: import convolve_cdef In [12]: timeit -n2 -r3 convolve_cdef.naive_convolve(f,g) 2 loops, best of 3: 1.56 s per loop
3. 効率的な配列参照
実はこの段階ではまだnumpy配列のアクセスがボトルネックになっています。これをCレベルにするには配列の型と配列の次元数を教えてやる必要があります。
def naive_convolve(np.ndarray[DTYPE_t, ndim=2] f, np.ndarray[DTYPE_t, ndim=2] g): … cdef np.ndarray[DTYPE_t, ndim=2] h = …
このチュートリアルではfとgに2次元のnumpy.int型を渡しているので、引数の型情報にDTYPE_t（defineしてあるので実際はnumpy.int_t）とndim=2を明示する必要があります。どうもバイト数が合っていればいけるみたいで、DTYPE_tの代わりにlongとしても動作しました。が、分かりにくくなるので普通にnumpy配列を宣言した時の型名に"_t"を付けたほうが良さそうです。
numpy.ndarrayに型情報を与えることによって、要素番号にint, unsigned intが使用されていればndarrayへのアクセスがCレベルに高速化されます。
実行結果
In [13]: import convolve_efficient In [14]: timeit -n2 -r3 convolve_efficient.naive_convolve(f,g) 2 loops, best of 3: 9.02 ms per loop
すげーー！！爆速！！今度はチュートリアル同様に圧倒的に速くなりましたね！（よかった...(^.^; ) 型情報をちゃんと与えることでCython+NumPyの真価が発揮されます。
最終的なコードはこちらです。といっても本家のチュートリアルのコードのまんまですが・・・
from __future__ import division import numpy as np cimport numpy as np DTYPE = np.int ctypedef np.int_t DTYPE_t def naive_convolve(np.ndarray[DTYPE_t, ndim=2] f, np.ndarray[DTYPE_t, ndim=2] g): if g.shape[0] % 2 != 1 or g.shape[1] % 2 != 1: raise ValueError("Only odd dimensions on filter supported") assert f.dtype == DTYPE and g.dtype == DTYPE cdef int vmax = f.shape[0] cdef int wmax = f.shape[1] cdef int smax = g.shape[0] cdef int tmax = g.shape[1] cdef int smid = smax // 2 cdef int tmid = tmax // 2 cdef int xmax = vmax + 2*smid cdef int ymax = wmax + 2*tmid cdef np.ndarray[DTYPE_t, ndim=2] h = np.zeros([xmax, ymax], dtype=DTYPE) cdef int x, y, s, t, v, w cdef int s_from, s_to, t_from, t_to cdef DTYPE_t value for x in xrange(xmax): for y in xrange(ymax): s_from = max(smid - x, -smid) s_to = min((xmax - x) - smid, smid + 1) t_from = max(tmid - y, -tmid) t_to = min((ymax - y) - tmid, tmid + 1) value = 0 for s in xrange(s_from, s_to): for t in xrange(t_from, t_to): v = x - smid + s w = y - tmid + t value += g[smid - s, tmid - t] * f[v, w] h[x, y] = value return h
実行時間比較
| time(s) | pythonとの相対倍率(倍) | |
| 元のPython | 2.13 | 1.0 |
| Cythonコード改変無し | 1.74 | 1.22 |
| 型情報追加 | 2.68 | 0.79 |
| ndarray以外だけに型情報追加 | 1.56 | 1.37 |
| ndarrayに詳細な型情報追加 | 0.00902 | 237 （！） |
チュートリアルではさらに、@cython.boundscheck(false)のデコレータを使用すると配列参照が配列の範囲内にあるかどうかのチェックを省略するのでさらに高速になると書いてあり、実行時間もさらに半分ぐらいになるようです。
が、私の環境ではチュートリアルのソースを丸コピしても何故か動かせませんでした・・・。
仕方がないので非常に簡単な関数を作って試してみたのですが効果があまり見られず、配列の参照をミスるとC言語で大変お世話になったsegmentation faultの文字と共にpythonを巻き込んで落ちました。ひぃぃ
ちなみに、boundscheck(false)のようなコンパイルオプションは他にもあって高速化に寄与するものもあるみたいです。興味のある方はこちらをどうぞhttp://wiki.cython.org/enhancements/compilerdirectives
自分のスキルではまだ下手にいじれなさそうです^^;
まとめ
- ループ変数は確実にcdef intしておく
- cimport numpyでCの領域でnumpyを扱う
- numpy配列もnp.ndarray[TYPE, ndim=N]の型情報を与える
- boundscheck(false)などのコンパイラオプションを与える
ここのチュートリアルでは登場しませんでしたが、他にも
- 関数をcdefにすることで呼び出しコストを低くする
- Pythonの関数の代わりにCの関数を呼び出す
- クラスをcdefにする
- prangeなどでGILを外して並列化
などの方法で高速化できるらしいです。覚えるのも大変そうですが夢が広がりますね！
それでは皆さんも素敵なCythonライフを！
今回使用したコードをgithubにアップしました
Pythonを高速化するCythonを使ってみた
突然ですが私はPythonが好きです。でもPythonは遅いです。
何が遅いかというと、致命的なことに四則演算が遅いです。でも他の動的型付け言語でスクリプト言語と呼ばれるPerl, Ruby, Javascript も C, Javaのようなコンパイルを行う静的型付け言語に比べれば圧倒的に遅いです（近年ではJavascriptのように著しく進歩した言語もあるので必ずしもそうだとは言えませんが）。
スクリプト言語が遅い原因の一つは、変数の型が指定されていないので型のチェックを毎回行う必要があるからです。この特性があるおかげ自動的に型を変換してオーバーフローを防いでくれるというメリットもあるのですが、どうしても静的型付け言語よりは速度を出すことができません。
ならPythonのコードに型指定を加えてコンパイルしちゃえばいいじゃん！というのがCythonです。正確にはPythonライクな文法で書いたコードをC/C++に変換してコンパイルします。噂では単純な計算だとPythonのコードを実行するより100倍以上も高速化することもあるらしい（！）ということで試してみました。
Cythonの基本的な使用方法は公式のドキュメントをご覧下さい
http://docs.cython.org/
日本語版
http://docscythonja.zouri.jp/
注意：以下のベンチマークはMacOS 10.7.3 MacBook Core2 Duo 2.26GHzでPythonとC/C++のコンパイラはMac標準のPython2.7.1, llvm-gccを使用しています。
C/C++はtimeコマンド、Python/Cythonはipythonのtimeitを利用して実行時間を測定しています。
なお、以下の記録はあくまで私の環境、私の実装での記録なので比較の結果は正しいというものではないです。記事の下の方に他の方のベンチマークへのリンクを載せましたのでそちらもご覧下さい。
単純な加算
まずは単純な加算を比較してみましょう。1~100000000までの数を足しあわせます。
C
$ gcc -O2 simple_add.c $ time ./a.out real 0m0.004s user 0m0.001s sys 0m0.003s
Python
In [1]: import test_py In [2]: timeit -n1 test_py.simple_add(100000000) 1 loops, best of 3: 8.39 s per loop
Python遅っ！Cは多分速すぎてちゃんと測定できてないです。
ではCythonを使ってみましょう。pyximportを使用するとCythonがimport時にコンパイルしてくれます。まずはPythonのコードを何も変更しないでただコンパイルしてみましょう．
Python/pyximport+.py
In [1]: import pyximport; pyximport.install(pyimport = True) In [2]: import test_py In [3]: timeit -n1 test_py.simple_add(100000000) 1 loops, best of 3: 4.27 s per loop
おお！何もしないで2倍近く速くなりました！これだけでも十分ありがたい結果です。
次はコードに何も手を付けずに拡張子をCythonの.pyxに変更してみます
Cython/pyximport+.pyx
In [1]: import pyximport; pyximport.install() In [2]: import test_cy In [3]: timeit -n1 test_cy.simple_add(100000000) 1 loops, best of 3: 4.36 s per loop
前の結果とほとんど変わらないですね。誤差の範囲でしょう。
さていよいよ本番です。pyxで型の指定を行います。cdefで型を指定することでCの変数を使うので速くなるそうです。
Cython/cdef
In [1]: import pyximport; pyximport.install() In [2]: import test_cy In [3]: timeit -n1 test_cy.simple_addc(100000000) 1 loops, best of 3: 1.19 us per loop
・・・( ﾟдﾟ)
Pythonと次元が違う、というか時間のオーダーが全然違います。Cの時間がtimeコマンドでは測定不能でしたが、これに迫る速さ？
if文と動的配列
次はif文を使用して10000000までの奇数を動的に配列に入れるという処理をしてみます。
C++は慣れていないですが、vectorを使って実装してみました。Pythonはlist.appendを使用します。
Python/list
In [1]: import test_py In [2]: timeit -n1 test_py.odd(10000000) 1 loops, best of 3: 2.24 s per loop
list.appendとかありえん！Python使いだったらリスト内包表記だろ常考！という声が聞こえてきますのでリスト内包表記でも試してみます。
Python/list内包表記
In [1]: import test_py In [2]: timeit -n1 test_py.odd_listcomp(10000000) 1 loops, best of 3: 1.64 s per loop
なるほど、内包表記にするとたしかに速くなります。しかし、さすがにC++と比べると全然といったところ。
Cythonでは問題なくリスト内表記が使用できるので、Pythonのコードに型指定だけ加えてコンパイルします。ここからはsetup.pyを用意してコンパイルを行います。
Cython/cdef+list内包表記
In [1]: import test_cy In [2]: timeit -n1 test_cy.odd_listcomp(10000000) 1 loops, best of 3: 288 ms per loop
C++に届きませんが、元のPythonよりかなり高速になりました。Cythonでもリスト内包表記が使えるのはありがたいです。
次はCythonでC++のvectorを呼び出してみます．実はCythonではC++の標準ライブラリが簡単に使えるようになっていて、deque, list, map, pair, queue, set, stack, vectorがサポートされているようです。vector
配列に対する計算
3つ目は、NumPyとも比較してみます。1~10000000までの数を格納した配列全てでlogを取ってから、逐次読み取ってsumを計算します。NumPy以外はlogとsumの計算を一度にやってしまってもいいのですが、NumPyに有利なようにlogとsumの行程を分けてあります。また、念のため計算結果を表示します。
まずはC++とPythonから。
C/vector+math.log
$ g++ -O2 logsum.cpp $ time ./a.out 1.51181e+08 real 0m0.742s user 0m0.515s sys 0m0.197s
Python/math.log
In [1]: import test_py In [2]: li = [i for i in xrange(1, 10000000)] In [3]: timeit -n1 test_py.logsum(li) 151180949.369 151180949.369 151180949.369 1 loops, best of 3: 3.48 s per loop
次はCythonでC++と同じようにvectorを使用します。また、Pythonのmathを使用しないでCのmath.hを使用する方が速いそうなのでそちらを利用します。
Cython/vector+cdef math.log
In [1]: import test_cy In [2]: li = [i for i in xrange(1, 10000000)] In [3]: timeit -n1 test_cy.logsum(li) 151180949.369 151180949.369 151180949.369 1 loops, best of 3: 850 ms per loop
Python/Numpy
In [1]: import test_py In [2]: li = [i for i in xrange(1, 10000000)] In [3]: timeit -n1 test_py.logsum_numpy(li) 151180949.369 151180949.369 151180949.369 1 loops, best of 3: 2.09 s per loop
うーん・・・単純にPythonで書くよりは高速ですが、Cython/vectorの結果を見てしまうと(^_^;)
NumPyはSIMD命令に対応しているらしいのでもっと速いと思っていたのにかなりガッカリです。
NumPyは慣れてしまえば行列に対する計算を分かりやすくコンパクトに書けるのですが、本気で速度を気にするならCythonでゴリゴリ書くべきなんでしょうかね・・・（本当はC++でゴリゴリ書くべき）
＊オマケ
CythonはNumPyにもサポートしているそうなのでどれだけ高速化するのか試してみます。
10000*10000の行列に0から順に入れて、全て+1した後に行ごとのsumで割って、行列全ての合計を求める。行ごとのsumで割ると行の合計は1.0になるはずなので、sum=1*1000になるはずです。
Python/Numpy
In [1]: import test_py In [2]: timeit -n1 test_py.dim2sum_numpy(1000) 1000.0 1000.0 1000.0 1 loops, best of 3: 869 ms per loop
Cython/Numpy+cdef int
In [1]: import test_cy In [2]: timeit -n1 test_cy.dim2sum_numpy(1000) 1000.0 1000.0 1000.0 1 loops, best of 3: 621 ms per loop
CythonはPython/NumPyのコードのループに使用する変数をint型に指定しただけです。劇的なスピードアップはありませんが、ほぼコードの改変無しにしては十分でしょう。
次はcimportでNumPyを使用します。こうするとNumPyの配列に高速にアクセスできるらしい。
Cython/cimport Numpy+cdef int
In [1]: import test_cy In [2]: timeit -n1 test_cy.dim2sum_pynumpy(1000) 1000.0 1000.0 1000.0 1 loops, best of 3: 630 ms per loop
あれ、変わらない？cimportで効果があるのは引数にNumPyの配列を受け取るときだけなのだろうか？
うーん、勉強不足ですいません(^_^;)
結果：
| time(s) | Pythonとの相対効率(倍) | |
|---|---|---|
| Python/NumPy | 0.869 | 1.0 |
| Cython/NumPy+cdef int | 0.621 | 1.4 |
| Cython/cimport NumPy+cdef int | 0.630 | 1.38 |
追記(3/10):
@lucidfrontier45さんに時間が変わらないのはnumpy.ndarrayをcdefしていないからだと教えて頂きました。
さらに、そもそもndarrayの配列の参照はfoo[i][j]よりもfoo[i,j]の方が速いことが分かったので、この2点を修正してやりなおしました。
結果（修正版）：
| time(s) | Pythonとの相対効率(倍) | |
|---|---|---|
| Python/NumPy | 0.350 | 1.0 |
| Cython/NumPy+cdef int | 0.222 | 1.58 |
| Cython/cimport NumPy+cdef int | 0.0431 | 8.12 |
Cythonを使わなくても2次元のndarrayの参照の方法を変えただけで速くなりました！これは今後も気を付けなければ・・・。
Cythonの方はndarrayをcdefしてやるとかなり速くなることが確認できました。今回のコードだとndarray.sum()を使用している部分が多いですが、配列の逐次参照部分が多ければもっと差がつくのではないかと思います。
アップしてあるコードも修正しました
他の方のベンチマーク：
High Performance Python tutorial v0.2 (from EuroPython 2011) 上の方のpdf
Pythonを高速化する様々な手法のベンチマークが載っている．
Python2.7 < PyPy1.5 < Cython < Cython+numpy < ShedSkin
の順番に高速だったらしい．ShedSkinはCythonと同様にPythonをC++に変換するものらしいけどまだ実験段階？
Speeding up Python (NumPy, Cython, and Weave)
Python < NumPy < Weave < Cython
WeaveというのはSciPyのモジュールでPythonのコードにインラインでCのコードを書けるものらしい．
結論
今回の結果と自分のスキル，可読性を考えると
Python << NumPy < Cython < Cython+NumPy
＊例外として互換性を完全に捨てる形でほぼCのCython
という感じだろうか．
NumPyとCythonを勉強していくうちに順位が変わる可能性は大いにありそうです
追記(2013/9/16):
かなり今更ですが、Cython公式のチュートリアルに沿って、NumPy+Cythonの組み合わせをさらに高速化する記事を書いていました。NumPyとCythonを組み合わせると爆速！
その記事の最後に触れていますが、Cythonのコンパイルオプションを変更すると、デフォルトの設定よりもさらに高速化させることができます。
残念ながらまだ和訳はされていないようですが、興味がある方は本家ドキュメントのこちらをどうぞ
http://wiki.cython.org/enhancements/compilerdirectives
SuffixArrayを作ってみた(Python)
前回の続き
自然言語処理はじめました - Ngramを数え上げまくるを見てSuffixArrayを使ったNgram取得をやってみたくなりました！
前回SuffixArray構築に必須の文字列のソートアルゴリズムのマルチキークイックソートが実装できたのでいよいよSuffixArrayを実装します。
SuffixArrayの構築や、これを使用した文字列の検索は単純な実装なら簡単です。自分はこんな感じに作りました
- SuffixArray構築
- 文字列内の各文字のポインタを保存（Pythonではポインタがないので要素番号）
- ポインタの位置から後ろの文字列でソートしてポインタを入れ替える
- 部分文字列の検索
- SuffixArray内から二分探索で検索文字列を検索
- 発見できたらその前や後ろにも検索文字列が存在しないかチェック
- 発見できた部分文字列の開始ポインタ位置を全て出力
- Ngram数え上げ
- SuffixArrayを端からngram見ていき、同じ文字列だったらカウント+。違う文字列になったら出力
#coding: utf-8 ''' 単純なSuffixArrayの実装 部分文字列の検索 Ngram数え上げ ''' def multikeysort(string, indices, cmpindex=0): """単純なマルチキークイックソート 文字列のインデックスを使用したソート""" if not indices: return [] pivot = indices[0] left=[] mid=[] right=[] #cmpindex番目の文字を比較 for index in (indices[i] for i in xrange(1, len(indices))): if string[index+cmpindex] < string[pivot+cmpindex]: left.append(index) elif string[index+cmpindex] == string[pivot+cmpindex]: mid.append(index) else: right.append(index) mid.append(pivot) #番兵"$"によるIndexError防止と計算量軽減 if not string[pivot+cmpindex] == '$' and not len(mid) < 2: mid = multikeysort(string, mid, cmpindex+1) left = multikeysort(string, left, cmpindex) right = multikeysort(string, right, cmpindex) return left + mid + right def get_suffixarray(string): """文字列からSuffixArrayを構築""" indices = xrange(len(string)) return multikeysort(string, indices) def binary_search(string, suffixarray, query): """部分文字列を見つけたときのmidを返す""" endindex = len(query) stringlen = len(string) min = 0 max = len(suffixarray)-1 while True: mid = (min+max) / 2 suffix = suffixarray[mid] #IndexError防止 if not suffix + endindex > stringlen: cmpstring = string[suffix:suffix+endindex] if max < min: return -1 elif cmpstring < query: min = mid + 1 elif cmpstring > query: max = mid - 1 #elif string[suffix][:endindex] == query: else: return mid def search_all(string, suffixarray, query): """部分文字列の位置を全て検索 見つかれば接尾辞要素番号""" endindex = len(query) mid = binary_search(string, suffixarray, query) left = mid right = mid + 1 while True: suffix = suffixarray[left] cmpstring = string[suffix:suffix+endindex] if cmpstring == query: yield suffix left -= 1 else: break while True: suffix = suffixarray[right] cmpstring = string[suffix:suffix+endindex] if cmpstring == query: yield suffix right += 1 else: break def count_ngram(string, suffixarray, n): """n-gramを数えるイテレータ""" #SuffixArrayからngram以上の文字列を取得するジェネレータ。"$"を除くので-1 cmp_strings = (string[i:i+n] for i in suffixarray if not i+n > len(string)-1) #初期化 count_string = cmp_strings.next() count = 1 for cmp_string in cmp_strings: if count_string == cmp_string: count += 1 else: yield (count_string, count) count_string = cmp_string count = 1 yield (count_string, count) #SuffixArrayの構築 >>> string = 'abracadabra$' >>> suffixarray = get_suffixarray(string) >>> for suffix in [string[i:] for i in suffixarray]: ... print suffix $ a$ abra$ abracadabra$ acadabra$ adabra$ bra$ bracadabra$ cadabra$ dabra$ ra$ racadabra$ >>> string = u"イカちゃんかわいいイカちゃん" + "$" >>> suffixarray = get_suffixarray(string) >>> for suffix in [string[i:] for i in suffixarray]: ... print suffix ... $ いいイカちゃん$ いイカちゃん$ かわいいイカちゃん$ ちゃん$ ちゃんかわいいイカちゃん$ ゃん$ ゃんかわいいイカちゃん$ わいいイカちゃん$ ん$ んかわいいイカちゃん$ イカちゃん$ イカちゃんかわいいイカちゃん$ カちゃん$ カちゃんかわいいイカちゃん$ #部分文字列の検索 >>> query = u'イカ' >>> print [i for i in search_all(string, suffixarray, query)] [9, 0] #Ngram数え上げ >>> for word, n in count_ngram(string, suffixarray, 2): ... print word, n ... いい 1 いイ 1 かわ 1 ちゃ 2 ゃん 2 わい 1 んか 1 イカ 2 カち 2
できたことにはできたけど、検索とngram数え上げの部分がかなり見苦しい・・・。番兵を上手く使えよって話なんですが、一文字ずつ比較すると遅くなりそうな気がしてしまってこんな感じにせざるを得なかったです。他にいい方法ないかな。
より実用的なものを目指すのであればSuffixArrayの圧縮や、効率的なソート、検索手法を組み合わせる必要があるかと思います。ちょっと調べたんですが、複雑すぎてサッパリ(-_-;)突き詰めるならデータ構造について知識がないとダメですね。
マルチキークイックソート
自然言語処理はじめました - Ngramを数え上げまくるを見てSuffixArrayを使ったNgram取得をやってみたくなりました！
SuffixArrayの構築には文字列のソートが必要なのでまずはここから始めます。
文字列のソートにはクイックソートを文字列のソートに応用したマルチキークイックソートを使います。アルゴリズムはこちらを参考にしました。
マルチキー・クイックソート(multikey-quicksort)で高速に文字列ソート - EchizenBlog-Zwei
まずは単純なクイックソートから
#coding: utf-8 def quicksort(seq): """単純なクイックソート 数値のソート""" if not seq: return [] pivot = seq.pop(0) left=[] right=[] for i in xrange(len(seq)): if seq[i] < pivot: left.append(seq[i]) else: right.append(seq[i]) left = quicksort(left) right = quicksort(right) return left + [pivot] + right >>> seq = [4,2,3,6,1,9,7,5,10] >>> print quicksort(seq) [1, 2, 3, 4, 5, 6, 7, 9, 10]
OK！リスト内包表記を使うとワンライナーでも書けるようですが分かりやすさ重視ということで^^;
次はこのクイックソートを参考にしながらマルチキークイックソート
>>> def multikeysort(seq, cmpindex=0): ... """単純なマルチキークイックソート ... 文字列のソート""" ... if not seq: ... return [] ... pivot = seq[0] ... left=[] ... mid=[] ... right=[] ... #cmpindex番目の文字を比較 ... for i in xrange(1, len(seq)): ... if seq[i][cmpindex] < pivot[cmpindex]: ... left.append(seq[i]) ... elif seq[i][cmpindex] == pivot[cmpindex]: ... mid.append(seq[i]) ... else: ... right.append(seq[i]) ... mid.append(pivot) ... #IndexError防止と計算量軽減 ... if not pivot[cmpindex] == '$' and not len(mid) < 2: ... mid = multikeysort(mid, cmpindex+1) ... left = multikeysort(left, cmpindex) ... right = multikeysort(right, cmpindex) ... return left + mid + right ... #英語 >>> string = "abracadabra$" >>> strings = [string[i:] for i in xrange(len(string))] >>> for suffix in strings: ... print suffix ... abracadabra$ bracadabra$ racadabra$ acadabra$ cadabra$ adabra$ dabra$ abra$ bra$ ra$ a$ $ >>> for suffix in multikeysort(strings): ... print suffix ... $ a$ abra$ abracadabra$ acadabra$ adabra$ bra$ bracadabra$ cadabra$ dabra$ ra$ racadabra$ #ユニコード >>> string = u"ソートは難しい$" >>> strings = [string[i:] for i in xrange(len(string))] >>> for suffix in strings: ... print suffix ... ソートは難しい$ ートは難しい$ トは難しい$ は難しい$ 難しい$ しい$ い$ $ >>> for suffix in multikeysort(strings): ... print suffix ... $ い$ しい$ は難しい$ ソートは難しい$ トは難しい$ ートは難しい$ 難しい$
とりあえず書けたというレベル^^;
今回は単純にleft, mid, rightの３つに分割しましたが、この処理はもっと効率の良い方法が存在するみたいです。
次はこのマルチキークイックソートを使ってSuffixArrayの構築、検索、Ngram取得の予定
