Fedora 23でchainerインストール


今,Fedora 23はgccのバージョンが5.3.1なので,そのままchainerをインストールしようとすると,コンパイルエラーが発生する.

原因は,以下の2点

  • CUDAがgcc 4.9までしか対応していないこと.
  • enumの扱いが変わったこと.

CUDAのgcc対応

コードを見ると,

/usr/local/cuda/include/host_config.h

の115行目に,GCCのバージョンが4.9より大きいとエラーになるような #error が書かれている.
仕方ないので,これをコメントアウトして次のようにするとコンパイルが通るようになる.

#error -- unsupported GNU version! gcc versions later than 4.9 are not supported!
   ↓
//#error -- unsupported GNU version! gcc versions later than 4.9 are not supported!

enumの扱い

元々enumはint型として扱えていたが,それができなくなっている.なので,

typedef enum enum_etype{
    a = 1
    b = 2,
    c = 3
} Etype;

int main(){
    Etype x = 0;
}

のようなコードはコンパイルエラーになる.これを防ぐにはgccのオプションに"-fpermissive" をつければ良い.

なので,chainerをインストールする際に,

alias g++="g++ -fpermissive"
pip install chainer

としたらいけた.正しいやり方かどうかわからないけど.

Pythonでbool値の引数をargparse


引数にTrue/Falseの値を渡してargparseでパースして使いたい.
ということでメモ.

ぐぐったらこういうやり方が出てきた.

def parse():
    parser = argparse.ArgumentParser(description='Argparseのテスト')
    parser.add_argument('--version', action='version', version='%(prog)s 1.0')

    parser.add_argument('--force', dest='force', action='store_true')
    parser.add_argument('--no-force', dest='force', action='store_false')
    parser.set_defaults(force=False)

    params = parser.parse_args()

    return vars(params)

この --force/--no-force の指定で force の値が True/False になる.デフォルトはFalse.

ただ,これだと

python hoge.py --force --no-force

みたいに両方指定することも可能なのでまずい.

そこで,mutually_exclusive_groupを使って,

def parse():
    parser = argparse.ArgumentParser(description='Argparseのテスト')
    parser.add_argument('--version', action='version', version='%(prog)s 1.0')

    group = parser.add_mutually_exclusive_group()
    group.add_argument('--force', action='store_true')
    group.add_argument('--no-force', action='store_false')
    parser.set_defaults(force=False)

    params = parser.parse_args()

    return vars(params)

とすることで,両方指定した時にエラーを出すことが出来る.

pep8準拠に


Python では, pep8 というコーディング規約(?)があって,
pep8 コマンドで,コードがそれに準拠しているかをチェックできる.

vim の syntastic と組み合わせると,ファイルを保存するたびに pep8 が実行され,
規約違反があった場合 vim の画面内の当該箇所に印がつく.
また,その行にカーソルを持って行くと,画面の一番下に,
どういう違反があるかが表示される.

例えば,変数名が短すぎる,関数名が hoge_hige() の形になっていない,コメント(docstring)がない,1行が長すぎるなど.

ひとまず自分の書いてるライブラリでエラーができるだけなくなるよう頑張る.

scikit-imageでの画像データの扱い


Pythonでの機械学習ライブラリとしては有名な scikit-learn があるが,
画像処理ライブラリは scikit-image がある.

最近Pythonでの画像処理の勉強をしていて気になったので,
scikit-image での画像の型の扱いについてメモ.

scikit-image ライブラリでは,画像の読み込みは,サブモジュールの skimage.data を使い,

import skimage.data

img = skimage.data.imread("lenna.jpg")

のように書く.

この時,読み込まれた画像は,uint8の型のarrayになっていて,
値は 0〜255 の間の整数値となる.

例えばグレースケール化したいと思った時,skimage.color.rgb2gray という関数を使うが,
これを行なうと,

import skimage.color
gimg = skimage.color.rgb2gray(img)

結果の gimg は,float64 型になっていて,
値は 0.0〜1.0 の浮動小数点型となる.

画像処理においては, 0.0〜1.0 の間の浮動小数点数として扱うのが,
メモリは食うけど便利.

この辺りの型の変換を行なうのが,

skimage.img_as_bool(img)
skimage.img_as_ubyte(img)
skimage.img_as_int(img)
skimage.img_as_uint(img)
skimage.img_as_float(img)

の関数.

それぞれ,
skimage.img_as_bool(img):
bool型 (True / False の2値)

skimage.img_as_ubyte(img):
uint8型(読み込んだ時の状態.[0,255]の範囲)

skimage.img_as_int(img):
int16型(uint8の[0,255]が[0,32767]にマップされる.負の値も格納可.)

skimage.img_as_uint(img):
uint16型([0,65535]の範囲)

skimage.img_as_float(img):
float64型(uint8の[0,255]が[0.0,1.0]にマップされる.負の値も格納可.)

となる.
なお,この範囲は, skimage.dtype_limits(img) によってタプルで取得可能.

3Dの素人がBundlerでVLFeatを使ってみた


この記事について

この記事はComputer Vision Advent Calendar 2012の12月10日分の記事として書いています.(一部はスタバにて,Mac Book Proで書いています.内容のBundler利用はLinux (Fedora 17) で行なっています.)

さて,今回の内容は,Bundlerを動かすこと,そしてその中で出てくるプログラムの一部にVLFeatというプログラムを利用する,という内容です.

Bundlerとは

Bundlerとは,Rubyのパッケージを管理するgemのラッパーで…ではなく,画像から3次元モデルを作るStructure from Motion (SfM)をするためのツールです.
オフィシャルのサイトにも

Bundler is a structure-from-motion (SfM) system for unordered image collections (for instance, images from the Internet) written in C and C++. An earlier version of this SfM system was used in the Photo Tourism project.

とあり,C/C++で実装されており,Photo Tourismというプロジェクトでも利用されていました.
で,なぜBundlerという名前かというと,このStructure from Motionの一連の処理の中でBundle Adjustmentという3次元と2次元の対応関係を計算する処理(ここがメイン?)があり,そこから名付けられているのだと思います.

基本的な処理の流れは,

  1. 特徴点を抽出する
  2. 複数の画像間で対応する特徴点を探索する
  3. 対応づいた特徴点をもとにBundle Adjustmentを実行する

となっています.
Bundlerのページからダウンロードできるコード(ver. 0.4)では,RunBundler.shというスクリプトにより,この一連の処理を行なっており,特にステップ1,2に対してはSIFT特徴を使って検出と対応関係の獲得を行なっています.

Bundlerの使い方自体は,満上師匠先生による解説が参考になります.
また,Bundle Adjustment自体の計算等については,岡谷先生のチュートリアル(もしくはコンピュータビジョン最先端ガイド)が参考になります.

オリジナルのBundlerの1,2のSIFTは,SIFT特徴を発明したD.Lowe先生の実装が利用されていますが,満上さんの話によるとR.Hess氏の実装(OpenSIFT)も簡単に利用できるとのことです.
今回,さらに別のSIFT実装であるVLFeatを利用してみようという話です.

ちなみに,Bundlerをコンパイルする時,新しいコンパイラでは基底クラスのコンストラクタ呼び出しの仕方でエラーが出ます.
src/BundlerApp.h の 620行目を消し,619行目に基底クラスのコンストラクタ呼び出しを追加してください.

VLFeatとは

VLFeatは様々な特徴量やアルゴリズムが実装されたライブラリです.

The VLFeat open source library implements popular computer vision algorithms including HOG,SIFTMSERk-meanshierarchical k-meansagglomerative information bottleneckSLIC superpixels, and quick shift. It is written in C for efficiency and compatibility, with interfaces in MATLAB for ease of use, and detailed documentation throughout. It supports WindowsMac OS X, and Linux. The latest version of VLFeat is 0.9.16.

HOGとかSIFTとかMSERとか色んな特徴量が実装されています.
(あれ…「OpenCVでええやん」という声が 聞こえてきた…)
まぁこれをダウンロードしてきてコンパイルすると,D.Lowe先生のSIFT実装の実行プログラムとほぼ同じ入出力のsiftという実行プログラムができます.

今回,このVLFeatで実装されているsiftプログラムをBundlerで使ってみます.

BundlerでVLFeatのsiftを使う

bundlerを展開したディレクトリ内のbinディレクトリにD.Lowe's siftの代わりにVLFeatのsiftを置くだけで済む…かと思いきや….

D.Lowe's siftとVLFeatのsiftの微妙な出力の違いでど嵌りしました.

だいぶ戦った結果,ステップ1で一旦出力するsift特徴量を記述したファイル(hoge.jpgが入力ならばhoge.key.gz)の形式が微妙に違うことがわかりました.

hoge.sift.gzをgunzipしてできたhoge.siftを見てみると,VLFeatに入っているsiftプログラムの出力は各行に

x y scale orientation 128次元がスペース区切り

が出力されていますが,D.Lowe's siftの出力には,先頭行に

特徴点数 次元数

の行があります.さらに,残りの各行は

20次元がスペース区切り
20次元がスペース区切り
20次元がスペース区切り
20次元がスペース区切り
20次元がスペース区切り
20次元がスペース区切り
8次元がスペース区切り

の形になっています.
Bundlerではこのファイルを開き,まず特徴点数を読み込んでから特徴点の配列を確保する仕様になっているので, VLFeatの入出力を少しいじる必要があります.
(Bundler側の読み込みコードを修正するという手もあるが,コードがあまりにもクソアレだったのでおすすめできない.)

sift入力の修正

D.Lowe's siftの場合,入力としては標準入力からpgmフォーマットの画像を読み込み,特徴量を標準出力へ出力する仕様になっているのに対し,VLFeatのsiftは,第1引数で与えられたhige/hoge.pgmを読み込み,hoge.keyをカレントディレクトリに出力する仕様になっている.

そのため,Bundlerのbin/ToSift.shファイルの36行目を次のように修正する.

修正前

echo "mogrify -format pgm $IMAGE_DIR/$d; $SIFT < $pgm_file > $key_file; rm $pgm_file; gzip -f $key_file"

修正後

echo sift_file=`echo $d | sed 's/jpg$/sift/'`
echo "mogrify -format pgm $IMAGE_DIR/$d; $SIFT $pgm_file; mv $sift_file $key_file; rm $pgm_file; gzip -f $key_file"

sift出力結果の修正

変換の修正ですが,特徴量ファイルはgz形式で圧縮されているので(ホントは上のToSift.shでgzipしなければ良いだけだけれど),

  1. gunzipで解凍
  2. 特徴量ファイルを修正(先頭行に特徴点数と次元数を追記)
  3. 各行を20次元ずつに分割
  4. gzipで圧縮

という手順になります.
まず,ステップ2,3の部分をpythonで書いてみると

#!/usr/bin/env python
import sys

i=0
buf=[]
for line in sys.stdin:
    i+=1
    l=line.split(" ")
    buf.append(" ".join(l[0:4]))
    l=l[4:]
    for j in range(6):
        buf.append(" ".join(l[j*20:(j+1)*20]))
    buf.append(" ".join(l[120:128]))

print i,"128"
print "\n".join(buf)

となります.
これをbundlerのbinディレクトリに,vlfeat2lowe.py として保存し,実行権限を与えます.

あとは,これに解凍/圧縮と,各ファイルに対して処理できるように ループを回して,

for img in `ls -1 $IMAGE_DIR | egrep ".jpg"`
do
        echo $img
        name=`basename $img .jpg`
        gunzip $IMAGE_DIR/${name}.key.gz
        $BASE_PATH/bin/vlfeat2lowe.py < $IMAGE_DIR/${name}.key > tmp.key
        mv tmp.key $IMAGE_DIR/${name}.key
        gzip -f $IMAGE_DIR/${name}.key
done

このコードをRunBundler.shの66行目あたりに突っ込めば動くようになるはずです.
ホントはもうちょっとスマートに書きたいが…ひとまず突貫工事ということで.
スマートに書きなおすには大工事が必要そうだったので.

結果の比較

D.Lowe's siftとVLFeatのsiftでBundlerを実行した結果の比較として,両手法で推定した3次元推定結果を以下に示す.
ただし,Bundlerの後処理として,密な3次元点群を得るために古川氏PMVS2を実行している.
可視化にはmeshlabを利用し,スクリーンショットを撮った.

まず入力画像はこちら

まずはD.Lowe's siftを用いた結果.
こっちのほうが圧倒的に処理は速かった.
結果の画像はこちら.

次にVLFeatのsiftを用いた結果.
検出されるキーポイント数が全然違う…閾値が違うのかな?
結果の画像はこちら.

違いが分からん…pmvs2するべきではなかった.
pmvs2前の頂点数を比較すると,D.Lowe's sift版では1459点に対してVLFeat版では1907点でした.
SIFTの検出キーポイント数が多い分,対応付けも多くとれたのかな.

終わりに

あーあ,完全にネタかぶっているブログを見つけてしまった….
zcatなんてコマンドがあったのか…

次やるときはもうちょっとスマートに書きたい.

scikit-learnを使ってPCAによる次元削減を7行で書いてみた


高次元のデータを解析する場合、データの分布などが見てもわからないので、PCAを使って次元削減をしてから可視化する場合があります。
PCAについてはこれまでにこのブログでも何度も触れているのでそちらを参照してください。

今回はフランスのINRIAがFundingしているscikit-learnというライブラリを使って7行で(空行含まず)書いてみました。言語はPythonです。

import numpy
import sklearn.decomposition

dim=2
data=numpy.loadtxt("data.csv", delimiter=",")
pca=sklearn.decomposition.PCA(dim)

result=pca.fit_transform(data);
numpy.savetxt("result.csv", result, delimiter=",")

たったこれだけで、data.csvというファイル(行ごとに多次元データがカンマ区切りで入っている)を読み込み、PCAで次元削減(2次元へ圧縮)したあとresult.csvというファイルに書き出しています。

白色化する場合は pca.whiten=True という行を、fit_transform の前に入れればOKだと思います。

PCAというより、numpy.loadtxt / numpy.savetxt 素晴らしい!という内容になりました。

matplotlibで3次元データのプロット


Pythonのグラフ描画用のライブラリであるmatplotlibを使って、3次元データのプロットを行ったときのメモ。

matplotlibで2次元のデータ( y = f(x) のようなもの )の散布図を描画するのは簡単で、matplotlibの中のpyplotを使って

import numpy
import matplotlib.pyplot as plt

xs=numpy.random.rand(1,100)
ys=numpy.random.rand(1,100)

plt.scatter(xs, ys)
plt.show()

でできる。結果はこのような感じ。

(もちろん先に関数fを定義しておく必要はある。上の例ではランダムな値を入れている)
要は、同じ数のデータ系列が2つあれば、それをx,y軸の値に取るような散布図を描画できる。

これを3次元でやりたい。要はx, y, zの3つの同じ数のデータが合った時の散布図を作りたい。
これはちょっと面倒で、mpl_toolkits.mplot3dというのを使う。

import numpy
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D

xs=numpy.random.rand(1,100)
ys=numpy.random.rand(1,100)
zs=numpy.random.rand(1,100)

fig=plt.figure()
ax=Axes3D(fig)

ax.scatter3D(xs, ys, zs)
plt.show()

これで下のような感じのプロットができる。

Tokyo.SciPy#4 出張版 (aka Kan.SciPy#1)で発表しました


※だいぶ前の話ですが

先日のTokyo.SciPy#4 出張版 (aka Kan.SciPy#1)で発表してきました。
Tokyo.SciPyはその名の通り、東京で行われているSciPy (Scientific Python)の勉強会です。
今回、特別に出張版として関西で開催されました。
そのため、名前はKan.SciPyとなっています。

僕は画像処理でPythonを(そしてNumPy及びSciPyを)使い始めていたので、 それに関して発表させていただきました。
資料はslideshareにアップしています。

発表するにあたって色々いじったので、だいぶPythonで画像処理をする事に関するノウハウが貯まったように思います。
それに関して、今後またこのブログで書いて行きたいと思います。

Pythonで主成分分析をするコード


2012/02/29 修正しました

以前C++からOpenCVで主成分分析するコードを書きましたが、今回はOpenCVも使わず、Pythonで実装してみました。

import numpy
import numpy.linalg
class PCA(object):
    """主成分分析
        eigenvectors:行ベクトルとして固有ベクトルが格納されている
        eigenvalues:固有値の列.大きい順にソートされている
        mean:平均ベクトル(行ベクトル)
    """
    def __init__(self, dim):
        self.dim=dim

    def calc(self,src):
        """主成分分析の計算を行う
        src:各データが行ベクトルとして格納されている
        """
        shape=src.shape

        # 中心化して転置する.以降,各データは列ベクトル
        mu = (sum(numpy.asarray(src,dtype=numpy.float)) / float(len(src)))
        src_m = (src - mu).T

        if shape[0] < shape[1]: # データ数 < 次元数のとき
            # PCA後の最大次元数の決定
            if shape[0] < self.dim: self.dim=shape[0]

            print>>sys.stderr,"covariance matrix ...",
            n=numpy.dot(src_m.T, src_m)/float(shape[0])
            print>>sys.stderr,"done"

            print>>sys.stderr,"eigen value decomposition ...",
            l,v=numpy.linalg.eig(n)
            idx= l.argsort()
            l = l[idx][::-1]
            v = v[:,idx][:,::-1]
            print>>sys.stderr,"done"

            vm=numpy.dot(src_m, v)
            for i in range(len(l)):
                if l[i]<=0:
                    v=v[:,:i]
                    l=l[:i]
                    if self.dim < i: self.dim=i
                    break
                vm[:,i]=vm[:,i]/numpy.sqrt(shape[0]*l[i])

        else: # データ数 >= 次元数のとき

            if shape[1] < self.dim: self.dim=shape[1]
            cov=numpy.dot(src_m, src_m.T)/float(shape[0])
            l,vm = numpy.linalg.eig(cov)
            idx= l.argsort()
            l = l[idx][::-1]
            vm = vm[:,idx][:,::-1]

        self.eigenvectors=vm[:,:self.dim]
        self.eigenvalues=l[:self.dim]
        self.mean=mu

    def project(self,vec):
        return numpy.dot(self.eigenvectors.T,(vec-self.mean).T).T

    def backproject(self,vec):
        return (numpy.dot(self.eigenvectors,vec.T)).T+self.mean

画像処理でPCAをする場合、データ数に対して次元数がとても大きい時が多いのですが、その場合、普通に共分散行列を作るとメモリが溢れてしまいます。 例:ベクトルが100000次元、データ数1000の場合、共分散行列の計算を普通にやると100000\times 100000の行列になってしまいます。 そこで、共分散行列 XX^T の代わりに X^TXの行列を計算するやり方があります。 例:このやり方だと、1000\times 1000の行列で済みます。 上のコードでは、入力したデータの次元数、データ数に応じてこの2つの方法を切り替えています。 使い方の例は次のようになります。

 def main():
    import numpy     # データの読み込み(各データは行ベクトルとして入っている)
    data = numpy.loadtxt("data.csv",delimiter=",")
    pca=PCA(2)     # PCAにより固有ベクトルなどを求める
    pca.calc(data)     # pca.eigenvectors, pca.eigenvalues, pca.meanで固有ベクトル、固有値、平均が得られる
    # 固有ベクトルは列ベクトルとして得られる
    # 固有空間への投影(行ベクトルとして結果が得られる)
    print>>sys.stderr,"projecting ...",
    p1=pca.project(data)
    print>>sys.stderr,"done"

    # 固有空間からの逆投影
    print>>sys.stderr,"back projecting ...",
    v1=pca.backproject(p1)
    print>>sys.stderr,"done"

if __name__=="__main__":
    main()

これもシリーズ物で、次回以降、確率的主成分分析やカーネル主成分分析について紹介して行きます。

[shogun 1] shogunを使ってみよう


最近Shogun - A Large Scale Machine Learning Toolbox という機械学習のツールボックスを使い始めたので,shogunについて紹介します.

まず紹介ページの英語を適当に日本語に訳してみます.

SHOGUNは大規模なカーネル法と特にサポートベクタマシン(Support Vector Machine: SVM)[1]に重点を置いている機械学習ツールボックスである.いろいろなSVM実装(OCAS [21], Liblinear [20], LibSVM [2], SVMLight, [3] SVMLin [4] and GPDT [5])に対するインターフェースになる汎用的なSVMオブジェクトを提供する.SVMは様々なカーネルと組み合わせることができる.このツールボクスは線形カーネル,多項式カーネル,ガウシアンカーネル,シグモイドカーネルだけでなく,Locality Improved [6], Fischer [7], TOP [8], Spectrum [9], Weighted Degree Kernel (with shifts) [10] [11] [12]なども提供する.また効率的なLINADD [12]最適化法も実装している.線形SVMに対して,COFFIN framework [22][23]は疎行列,密行列他様々なデータタイプの混合に対してもオンデマンドで処理ができる.さらにSHOGUNは様々なpre-computedなカーネルで動くようになっている.重要なものとしては,サブカーネルの重み付き線形和による結合カーネルがある.この各サブカーネルは必ずしも同じドメインである必要はない.最適なサブカーネルの重みはMultiple Kernel Learning (MKL) [13][14][18][19]で学習できる.現在SVMは,1クラス,2クラス,多クラス分類,回帰問題を扱える.SHOGUNは線形判別分析(Linear Discriminant Analysis: LDA),線形計画マシン(Linear Programming Machine: LPM),パーセプトロン((Kernel) Perceptrons),隠れマルコフモデル(Hidden Markov Model)も実装されている.入力データは様々な型の疎・密行列を利用することができる.平均を引くなどの一連の前処理を各入力データに適用することもできる.
SHOGUNはC++で実装されており,Matlab, R, Octave, Python, Java, Ruby向けのインタフェースが用意されており,機械学習のオープンソースソフトウェアとして公開されている.

長々と書きましたが,重要なのは以下のポイントになります.

  • 様々なSVM実装が利用できる
  • 様々なカーネルを使える
  • SVM以外にもLDAとかHMMとかいろいろ使える
  • C++/Matlab/R/Octave/Python/Ruby/Javaなどで使える
  • オープンソースソフトウェアである

関西CV・PRML勉強会で紹介した時のスライドはこちらにあります.

僕はこのツールボックスを,Pythonから利用しはじめました.スライド中ではapt-getやMacPortsでインストールできると言っていますが,サンプルコードなどを入手したい場合は,gitを使ってソースコードをダウンロードし,コンパイルしてインストールするのがお勧めです.

gitがインストールされているならば

git clone git://github.com/shogun-toolbox/shogun.git
git clone git://github.com/shogun-toolbox/shogun-data.git

のコマンドで,shogunディレクトリとshogun-dataディレクトリができます.shogunディレクトリにはshogunのソースコードとサンプルコードが,shogun-dataディレクトリにはサンプルのデータが入っています.

ソースコードをダウンロードしたら,

cd shogun/src
./configure
make

し,管理者権限で

make install

することでインストールすることができます.
./configureの前に,swigがインストールされていないとpythonなどからの利用ができないので,先にインストールしておく必要があります.これはapt-getやMacPortsなどでインストールすることができます.

サンプルコードを動かすためには,サンプルデータのディレクトリの設定をする必要があります.shogunディレクトリ中で

rmdir data
ln -s ../shogun-data data

としてリンクを張っておく,もしくはshogun-dataの中身をshogun/dataへ移動すると完了です.

次のエントリではPython版のサンプルコードを動かすためのノウハウについて書きます.