Octaveでdeformable part modelを使ってみる (2)

著者のページから最新のコードをダウンロードして展開します。
ディレクトリ内でoctaveを立ち上げ、まずmexでCファイルをコンパイルします。

$ cd voc-release4
$ octave
octave:1> mex dt.cc
octave:2> mex features.cc
octave:3> mex fconv.cc
octave:4> mex getdetections.cc
octave:5> mex resize.cc
resize.cc: In function ‘void resize1dtran(double*, int, double*, int, int, int)’:
resize.cc:46: error: ‘sy’ was not declared in this scope

resize.ccのコンパイルでエラーがでます。これは多分'sy1'の誤りですが、いずれにせよただのassertなのでコメントアウトしてもいいでしょう。

いよいよ、サンプルスクリプトのdemo.mを走らせてみます。バスの画像が表示されました。

しかし、次へ進むとエラーが出て落ちてしまいました。
HOGpicture.mの12行目で読んでるimrotateが原因とのこと。どうやらoctaveのimrotateでは補間方法をちゃんと書かないといけないようです。
以下のように修正します。

12   bim(:,:,i) = imrotate(bim1, -(i-1)*20, 'bilinear', 'crop');

すると、モデルが表示されるところまで行けるようになりました。論文とかでよく見るやつですね。

先へ進むとまた落ちてしまいます。今度はreduceboxes.mの33行目で代入の次元がおかしいといわれています。ただ、これはそもそもIの値が空(boxが見つからない)時に処理が進んでいるのが原因です。
(変数が空の時の処理はoctavematlabで異なるのでしょうか?matlab環境が手元にないので分かりません...)
以下のように、空の時は処理を飛ばすようにします。

 18   % process boxes for component i
 19   I = find(boxes(:,end-1) == i);
 20   if isempty(I)
 21     continue
 22   end
 23   tmp = boxes(I,:);

これでようやく検出結果が表示されるようになりました。

Bounding boxはこんな感じ。

しかし安心したのも束の間、次のpersonの例のbox描画部でまた落ちてしまいます。
showboxes.mの76行目のlineで怒られています。
octaveのlineでは、座標が配列でまとめて与えられるとまずいようです。なので、今回のように検出対象が複数だと問題が起こってしまいます。
仕方ないので地道にfor文で書き換えます。

for j = 1:size(x1,1)
  line([x1(j) x1(j) x2(j) x2(j) x1(j)], [y1(j) y2(j) y2(j) y1(j) y1(j)], 'color', c, 'linewidth', cwidth, 'linestyle', s);
end

これでめでたく3例ともきちんと動くようになりました。(まだ表面化していない不具合はあるかもしれませんが...)

追記

fconvはマルチスレッドのものをコンパイルしたほうが速いです。

octave:1> mex fconvMT.cc -o fconv

Octaveでdeformable part modelを使ってみる (1)

Deformable part model*1は近年特に注目を浴びている一般物体検出手法で、デファクトスタンダードとしての地位を固めつつあります。人気の理由として、認識精度のよさや理論的な面白さなどがあげられますが、やはり開発者がしっかりしたソースコードを提供していることが大きいと思います。
http://www.cs.brown.edu/~pff/latent/

ただ、このコードはmatlab用のものです。一般人は普通matlabなんか持っていません。。
まあ最近ではpythonの実装*2などもあるようですが、やはりできればまず開発者のコードを動かしてみたいですよね。

なので、フリーのmatlabクローンであるoctaveでサンプルを動かしてみます。
Octaveはクローンといいつつmatlabとの互換性はかなり怪しかったのですが、最近はだいぶマシになってきたようです。ただ、まだ微妙に仕様が違うところが多いので注意が必要です。

Octaveのインストール

http://www.gnu.org/software/octave/からoctaveソースコードを落としてきます。3.4系以降ならいいと思います。私は3.4.0を使いました。
ソースコードを展開してconfigureします。オプションはいろいろありますが、画像を扱うためにはmagickのライブラリを指定する必要があります。ImageMagickを使う場合はこんな感じ。(先にImageMagickgnuplotは入れておいてください)

$ cd octave-3.4.0
$ ./configure --with-magick=ImageMagick

他のオプションはお好みで。(configureスクリプトの中にいろいろ説明があります)
多分いろいろライブラリが足りないと怒られるので、適宜入れてください。大抵のものはパッケージ管理ツールで入りますが、パフォーマンスを上げるためにblaslapackあたりは自分でソースからコンパイルしてもいいかもしれません。

configureが通ったらmake、make installでOK。make中になんか足らないと怒られることもありますが、その時も適宜インストールしてconfigureからやりなおしてください。

imageパッケージのインストール

Octaveに、画像処理のパッケージを追加で入れます。octave-forgeからimageパッケージをダウンロードします。ダウンロードした場所でoctaveを起動し、pkg installすると入ります。

$ octave
octave:1> pkg install image-1.0.15.tar.gz

これでoctaveの準備は整いました。

convertと標準入出力

そういえば、ImageMagickのconvertで標準入出力ってどうやって使うんだろうと気になったのでやってみました。といっても至極簡単で、

$ convert -strip XXX.jpg ppm:- | ./compute_gist

のようにハイフンで指定するだけのようです。明示的にフォーマットを指定したい時はその前に拡張子を記述します。これでppm以外のファイルも一応compute_gistで特徴抽出ができます。
(まあ、この場合OpenCVなり使ってcompute_gistのソースをいじった方がいいとは思いますが)

標準入力も同じように使えるので、htmlから動的に呼び出したりする場合には便利そうです。

画像サイズと認識率

GISTは比較的計算コストが小さい特徴量ですが、やはりそれなりに時間がかかります。手っ取り早く処理を軽くするには画像サイズを小さくすることが考えられますが、認識精度との兼ね合いが気になるので調べてみます。

ベンチマークは、前回と同じMITのシーン画像データセットを用います。1クラスあたり100サンプルの学習画像をランダムに選び、認識率を算出します。これを10回繰り返し平均をとります。識別器はSVMを用いました。
結果は以下のとおりです。

画像サイズ 32x32 64x64 128x128 192x192 256x256
認識率(%) 76.1 79.8 81.5 81.0 80.4

ということで、128x128が最もよい結果となりました。一般的なサムネイル画像くらいの大きさがあればOKなようです。もちろんデータの性質にもよりますが、あまり大きすぎると逆に良くないようです。

GISTで画像認識

せっかくなのでGISTで画像認識をやってみます。
GISTは画像全体のシーン認識を目的に設計されており、このタスクに特に適しているとされています。

まず、GISTの開発者が提供している、8クラスのシーン画像データセットで試してみます。
以下のページからImages.zipを落としてきて展開します。
http://people.csail.mit.edu/torralba/code/spatialenvelope/

一つのディレクトリの中に全ての画像が置かれていますが、扱いやすいようにクラスごとに分けることにします。

$ cd spatial_envelope_256x256_static_8outdoorcategories
$ find . -name "*jpg" | xargs -l -i basename {} | cut -d'_' -f1 | sort | uniq | xargs mkdir
$ mv mountain*jpg mountain
$ mv opencountry*jpg opencountry
$ mv forest*jpg forest
$ mv coast*jpg coast
$ mv tallbuilding*jpg tallbuilding
$ mv street*jpg street
$ mv insidecity*jpg insidecity
$ mv highway*jpg highway

例のコードはppm/pgmしか読めないので、ImageMagickで変換します。
ヘッダにコメント領域が残ってるとおかしくなることがあるので-stripをつけます。このデータセットの画像は256x256のサイズにもともと揃えてあるので、サイズに関してはこのままで大丈夫です*1

$ find . -name "*jpg" | sed 's/\.jpg//g' | xargs -l -i convert -strip "{}.jpg" "{}.ppm"
$ find . -name "*jpg" | xargs rm

準備ができたので、データセットの特徴を抽出します。perlで適当にスクリプトを書いてみました。

#!/usr/local/perl

use strict;
use List::Util;

my $Ntrain=100;
my $imgdir="/home/XXX/spatial_envelope_256x256_static_8outdoorcategories";
my $cmd="/home/XXX/lear_gist-1.1/compute_gist";

opendir(IMGDIR, $imgdir) or die($imgdir);
my @classes = readdir(IMGDIR);
close(IMGDIR);

open(TRAIN, ">train.txt");
open(TEST, ">test.txt");

foreach my $class(@classes){
    next if ($class =~ /^\./ || -f $class);

    opendir(CLASS, "$imgdir/$class");
    my @files = readdir(CLASS);
    closedir(CLASS);

    @files = List::Util::shuffle @files;

    my $count = 0;
    foreach my $file(@files){
        next if( -d "$imgdir/$class/$file");
        print "$imgdir/$class/$file\n";

        open(GIST, "$cmd $imgdir/$class/$file |") or die("$cmd");
        my $gist = <GIST>;
        close(GIST);
        next if(!$gist);

        if($count<$Ntrain){
            print TRAIN "$file $class $gist";
        }else{
            print TEST "$file $class $gist";
        }
	$count++;
    }
}

close(TRAIN);
close(TEST);

$imgdir以下におかれたディレクトリから、各カテゴリの学習サンプルを一定数サンプリングし、残りをテストサンプルとしています。
なんかオーバーヘッドで時間とられているような気もしますがとりあえず無視。

無事特徴もとれたので、RのSVMで識別をやってみます。
kernlabパッケージが必要なのでインストールします。

> install.packages('kernlab')

作った特徴ファイルを読み込んでSVMを実行します。

> library('kernlab')
> train <- read.table('train.txt')
> test <- read.table('test.txt')
> ot8.svm <- ksvm(train[,2]~.,data=train[,3:962],cross=5)
> ot8.pre <- predict(ot8.svm,test[3:962])
> (ot8.tab <- table(test[,2],ot8.pre))
> mean(diag(ot8.tab)/rowSums(ot8.tab))

[1] 0.8074061

という感じで、識別率は約80%となりました。本当はサンプルを入れ替えながら何回か試行するべきですが。

ついでにRandomForestも試してみました。パラメータはデフォルトで。

> install.packages('randomForest')
> library('randomForest')
> set.seed(20)
> ot8.rf <- randomForest(train[,2]~.,data=train[,3:962])
> ot8.pre <- predict(ot8.rf,test[3:962])
> (ot8.tab <- table(test[,2],ot8.pre))
> mean(diag(ot8.tab)/rowSums(ot8.tab))

[1] 0.7760976

ちなみに、前回やった判別分析だと

library('MASS')
set.seed(20)
ot8.lda <- lda(train[,2]~.,data=train[,3:962])
ot8.pre <- predict(ot8.lda,test[3:962])
(ot8.tab <- table(test[,2],ot8.pre$class))
mean(diag(ot8.tab)/rowSums(ot8.tab))

[1] 0.4530859

けっこう次元が大きいので、正則化項を入れてやらないと無理っぽいです。

*1:サイズも変えたいときは convert -strip -geometry 256x256! "{}.jpg" "{}.ppm" のようにしてやるとよいでしょう。!がないと長辺を256ピクセルにしてしまうので注意

GIST特徴抽出

GISTは代表的な大域的画像特徴量(画像全体から求まる特徴ベクトル)の一つです。コンピュータビジョン、特に一般物体認識の大家であるA.Torralbaらによって開発された特徴量*1で、現在でもしばしば用いられます。コピー画像検出タスクなどでは、bag-of-words (BoW)に匹敵する性能が得られることも報告されています*2。画像特徴とってみたいけど、BoWはなんかめんどくさそう・・・という人にはおすすめです。

本家のページでMatlabのコードが提供されています。
http://people.csail.mit.edu/torralba/code/spatialenvelope/

また、これまた有名なINRIAの研究グループがCのコードを提供しています。
http://lear.inrialpes.fr/software (一番下)
こちらを使ってみることにします。

コンパイルにはfftw3が必要。READMEに従ってhttp://www.fftw.org/からソースを落とす。普通にconfigure, make, make installでよいはず。以下のオプションは少なくとも必要。

$ cd fftw-3.3
$ ./configure --enable-float --enable-shared
$ make
$ make install

無事インストールできたら、GISTソースコードディレクトリへいってmake。

$ cd lear_gist-1.1
$ make

ディレクトリにあるサンプル画像の特徴をとってみます。

$ ./compute_gist ar.ppm
0.0579 0.1926 0.0933 ・・・ 0.0563 0.0575 0.0640

デフォルトだと、960次元のベクトルになっているはずです。画像全体を4x4のグリッドに分割し、各小領域で20方向のフィルタバンクの反応を計算します。これを各色について行うので、合計で4x4x20x3=960次元となります。パラメータはコマンドライン引数で変えられるようですが、特に理由がなければデフォルトでよいと思います。

なお、READMEに記述されているように、特徴抽出には注意が必要です。
===
The code works on non-square images, but it is of little interest:
descriptors computed on images of different sizes are not comparable.
===
つまり、特徴抽出を行う画像はあらかじめ同じサイズに揃えておく必要があります。
一般的には、適当な大きさの正方形に揃えるようです。(256x256とか)

*1:Aude Oliva, Antonio Torralba, Modeling the shape of the scene: a holistic representation of the spatial envelope, IJCV, Vol.42(3), pp.145-175, 2001.

*2: Matthijs Douze, Hervé Jégou, Harsimrat Sandhawalia, Laurent Amsaleg, Cordelia Schmid, Evaluation of GIST descriptors for web-scale image search, CIVR'09, 2009

CentOSにRを入れてみた

ここから必要なrpmパッケージを落とす。
まずR本体のインストール。

$ rpm -ivh R-2.13.0-2.el6.rf.x86_64.rpm
エラー: 依存性の欠如:
	libtcl8.5.so()(64bit) は R-2.13.0-2.el6.rf.x86_64 に必要とされています
	libtk8.5.so()(64bit) は R-2.13.0-2.el6.rf.x86_64 に必要とされています

必要なものをいろいろ言われるので適宜yumでインストールする。意外とめんどくさい。

$ yum install tcl tk

R-devel, libRmath, libRmath-develも入れる。

せっかくなので判別分析でもやってみる。かの有名なフィッシャーのあやめデータセットを使う。これはデフォルトで入っている。

> iris
    Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
1            5.1         3.5          1.4         0.2     setosa
2            4.9         3.0          1.4         0.2     setosa
3            4.7         3.2          1.3         0.2     setosa
4            4.6         3.1          1.5         0.2     setosa
5            5.0         3.6          1.4         0.2     setosa
~
145          6.7         3.3          5.7         2.5  virginica
146          6.7         3.0          5.2         2.3  virginica
147          6.3         2.5          5.0         1.9  virginica
148          6.5         3.0          5.2         2.0  virginica
149          6.2         3.4          5.4         2.3  virginica
150          5.9         3.0          5.1         1.8  virginica

3クラス、50サンプルずつデータが並んでいる。
LDAを使うにはMASSパッケージが必要。

> library(MASS)
> (Z.lda<-lda(iris[,5]~.,data=iris[,1:4]))
Call:
lda(iris[, 5] ~ ., data = iris[, 1:4])

Prior probabilities of groups:
    setosa versicolor  virginica 
 0.3333333  0.3333333  0.3333333 

Group means:
           Sepal.Length Sepal.Width Petal.Length Petal.Width
setosa            5.006       3.428        1.462       0.246
versicolor        5.936       2.770        4.260       1.326
virginica         6.588       2.974        5.552       2.026

Coefficients of linear discriminants:
                    LD1         LD2
Sepal.Length  0.8293776 -0.02410215
Sepal.Width   1.5344731 -2.16452123
Petal.Length -2.2012117  0.93192121
Petal.Width  -2.8104603 -2.83918785

Proportion of trace:
   LD1    LD2 
0.9912 0.0088 

識別実験もやってみる。奇数行を学習データ、偶数行をテストデータとする。見やすいようにクラスラベルも変える。

> iris.train<-iris1[even.n,]
> iris.lab<-c(rep("S",50),rep("C",50),rep("V",50))
> iris1<-data.frame(iris[,1:4],Species=iris.lab)
> even.n<-2*(1:75)-1
> iris.train<-iris1[even.n,]
> iris.test<-iris1[-even.n,]
> Z.lda<-lda(Species~.,data=iris.train)

学習データの識別結果はpredictで見れる。$classはクラス、$posteriorは事後確率、$xは判別特徴のスコア。
識別結果はtableで。

> predict(Z.lda)$class
 [1] S S S S S S S S S S S S S S S S S S S S S S S S S C C C C C C C C C C V C C
[39] C C C C C C C C C C C C V V V V V V V V V V V V V V V V V C V V V V V V V
Levels: C S V
> predict(Z.lda)$posterior
               C            S            V
1   2.063321e-21 1.000000e+00 6.123882e-42
3   1.241843e-18 1.000000e+00 2.385305e-38
5   5.002164e-22 1.000000e+00 1.218263e-42
7   4.894577e-18 1.000000e+00 5.109531e-37
9   1.138098e-14 1.000000e+00 3.454692e-33
~
> predict(Z.lda)$x
           LD1         LD2
1    7.9226230 -0.21852960
3    7.2987565  0.31563262
5    8.0475972 -0.46728206
7    7.0862313 -0.35174617
9    6.4034581  1.01764732
~
> table(iris.train[,5],predict(Z.lda)$class)
   
     C  S  V
  C 24  0  1
  S  0 25  0
  V  1  0 24

未知データを識別するには、以下のようにpredictを使う。

> Y<-predict(Z.lda,iris.test[,-5]) #使うのはテストデータの第5列以外
> table(iris.test[,5],Y$class) 
   
     C  S  V
  C 24  0  1
  S  0 25  0
  V  2  0 23