人工知能に関する断想録

■[目標] 2010年目標 20:49

• 人工知能に関する三テーマの深耕

特に興味を持っている人工知能に関する三テーマ（2009/3/23）をさらに深耕していきたい。いろんな題材はそろっているのだけど、どこからどう手をつけてどう形にしていくかがけっこう難しい。発明家を目指すならやっぱりモノを作らないと！今年は、去年の最後にやった3日で作る高速特定物体認識システム（2009/10/18）のような形式でまとめていこうかなと考え中。テーマに関連した技術や論文の紹介などは他の人にも役立つかもね。月に1つは記事にまとめたいとこだけど・・・お勤め疲れでなかなか研究時間がとれない。身体も鍛えないと。

• パターン認識と機械学習の上下を読破

昨年、「パターン認識と機械学習」を購入した（2009/5/7）のだが、すごく難しくて進捗があまり捗々しくなかった。関連する数学の修得も含めて、今年こそ読破するぞ。さっき、下巻もAmazonで購入した。8000円もしたので無駄にはできないぞ、これは！

• TOEICで900点

TOEICは2年前に受けてそれっきりだったけど、今年は3月、7月、11月に受験して英語力がどれくらい上がったか客観的に確かめようと思う。「英語上達完全マップを10ヶ月やってみた」の記事で知った英語上達完全マップにしたがって昨年末から訓練を初めている。音読パッケージと瞬間英作文をやってるのだけど、自分に欠けてたのはこれだ！と目から鱗だった。英語の修得は、海外でも困らないようにってのもあるけど、洋ゲー（最近はRisen, Dragon Age: Originsが面白い）をストレスなくプレイしたいってのもある（笑）これはモチベーション高いぞー。

その他の技術修得とか実験したいことはいろいろあるけど、主要目標はこの3つかな。昨年は、心の燈火（2005/6/13）がかき消えてしまう危機に見舞われた。実際、それに備えて心理学を学んでたのに・・・あまり役にたたなかった（むしろたったのか？）。今年は心機一転、モチベーションを維持してがんばるぞ！

■[目標] 2009年総括 22:26

2009年目標（2009/1/1）のつづき。

今年の目標はほとんど未達になってしまった・・・かろうじて画像処理・画像認識の技術を少し深耕したくらいかなぁ。もう、iPhoneアプリの作成もゲーム作成もずいぶん長い間ご無沙汰になってしまった。来年からは別の計画があるのでゲーム作成はもうやらないかもしれない。

今年は、死を身近に感じた非常につらい一年だった。来年はよい年にしたい。

O God, give us serenity to accept what cannot be changed,

courage to change what should be changed,

and wisdom to distinguish the one from the other.

Reinhold Niebuhr

他人と過去は変えられない、変えられるのは自分と未来だけだ。

■[OpenCV] 3日で作る高速特定物体認識システム (7) 最近傍探索の高速化 22:21

3日で作る高速特定物体認識システム (6) 線形探索を用いた特定物体認識（2009/11/22）のつづきです。今回がこのシリーズの最終回です。

前回の線形探索は遅すぎるので最近傍探索を高速化します。これで表題の高速特定物体認識システムができあがります。高速化にはいくつかの方法がありますが、物体モデルデータベースをなんらかのデータ構造にあらかじめ格納しておくというのがポイントです。今回は、資料でも述べられているkd-treeとLocality Sensitive Hashing (LSH)という手法を試してみます。kd-treeは木構造、LSHはハッシュでデータを構造化（インデキシング）します。kd-treeは、厳密な最近傍を求めますが、LSHは近似最近傍検索と呼ばれ、厳密な最近傍は求められない代わりに計算を大幅に高速化できます。

資料では、ANN (Approximate Nearest Neighbor）というライブラリのkd-treeを使っています。ANNのkd-treeは近似最近傍探索も行えるように拡張しているようです。OpenCVには、kd-tree（厳密版）とLSHは実装されています。OpenCVのkd-treeとLSHは、関数APIが似た形で整備されているので使うのは簡単です。

kd-tree

kdtree_recognition.cpp

まずはkd-treeです。えー、kd-treeとは、・・・というのは面倒なのでググってください。高次元データを効率よく検索するための木構造です。一般的には数十次元くらいまで高速に検索できると言われていますが、それ以上になると次元の呪いにより検索効率は落ちてきます。今回は、局所特徴量128次元なのですが線形探索に比べて大幅に高速化できました。

    // キーポイントの特徴ベクトルをobjMat行列にロード
    cout << "物体モデルデータベースをロードします ... " << flush;
    vector<int> labels;     // キーポイントのラベル（objMatに対応）
    vector<int> laplacians;  // キーポイントのラプラシアン
    CvMat* objMat;           // 各行が物体のキーポイントの特徴ベクトル
    if (!loadDescription(DESC_FILE, labels, laplacians, objMat)) {
        cerr << "cannot load description file" << endl;
        return 1;
    }
    cout << "OK" << endl;

    // 物体モデルデータベースをインデキシング
    cout << "物体モデルデータベースをインデキシングします ... " << flush;
    CvFeatureTree* ft = cvCreateKDTree(objMat);  // objMatはコピーされないので解放してはダメ
    cout << "OK" << endl;

kd-treeを作るにはcvCreateKDTree()関数を使います。引数は、1行に1データの行列です。物体モデルデータベースの全特徴量をobjMat（全特徴量数x128次元の行列）に展開しています。objMatのデータはコピーされずポインタをkd-treeに格納しているのでobjMatは最後まで解放してはダメです。

    // クエリのキーポイントの特徴ベクトルをCvMatに展開
    CvMat* queryMat = cvCreateMat(queryDescriptors->total, DIM, CV_32FC1);
    CvSeqReader reader;
    float* ptr = queryMat->data.fl;
    cvStartReadSeq(queryDescriptors, &reader);
    for (int i = 0; i < queryDescriptors->total; i++) {
        float* descriptor = (float*)reader.ptr;
        CV_NEXT_SEQ_ELEM(reader.seq->elem_size, reader);
        memcpy(ptr, descriptor, DIM * sizeof(float));  // DIM次元の特徴ベクトルをコピー
        ptr += DIM;
    }

    // kd-treeで1-NNのキーポイントインデックスを検索
    int k = 1;  // k-NNのk
    CvMat* indices = cvCreateMat(queryKeypoints->total, k, CV_32SC1);   // 1-NNのインデックス
    CvMat* dists = cvCreateMat(queryKeypoints->total, k, CV_64FC1);     // その距離
    cvFindFeatures(ft, queryMat, indices, dists, k, 250);

クエリの特徴ベクトル群を行列（CvMat）に展開して各特徴ベクトルの最近傍をkd-treeから検索する部分です。cvFindFeatures()という関数を使います。第1引数は、kd-treeです。第2引数はクエリ行列（各行がそれぞれクエリ画像の局所特徴量のベクトル）です。クエリ行列を与えると各特徴ベクトルのk-NNをまとめて計算できるようになってます。クエリが1つの特徴ベクトルの場合でもいったんCvMatに展開しないとダメなようですね。indicesとdistsは、クエリ行列と同じ行数を持つ行列で、クエリの各特徴ベクトルのk-NNのインデックスとクエリからの距離がそれぞれ格納されます。今回はk=1で1-NN（最近傍）の探索なので列数は1の縦に長い行列になります。k-NNを探索したい場合は、k列になります。小さなデータで試してみたところ距離はユークリッド距離で計算されているようです。

LSH

lsh_recognition.cpp

次にLSHを試してみます。実はLSHの実装はOpenCV2.0から追加されたようでマニュアルにも説明が書かれていませんでした。そんなわけでちょっと時間をかけてOpenCV2.0/src/cv/cvlsh.cppを解読して使い方を調べてみました。

    // キーポイントの特徴ベクトルをobjMat行列にロード
    cout << "物体モデルデータベースをロードします ... " << flush;
    vector<int> labels;      // キーポイントのラベル（objMatに対応）
    vector<int> laplacians;  // キーポイントのラプラシアン
    CvMat* objMat;           // 各行が物体のキーポイントの特徴ベクトル
    if (!loadDescription(DESC_FILE, labels, laplacians, objMat)) {
        cerr << "cannot load description file" << endl;
        return 1;
    }
    cout << "OK" << endl;

    // 物体モデルデータベースをインデキシング
    cout << "物体モデルデータベースをインデキシングします ... " << flush;
    CvLSH* lsh = cvCreateMemoryLSH(DIM, 100000, 5, 64, CV_32FC1);
    cvLSHAdd(lsh, objMat);
    cout << "OK" << endl;
    cout << "LSH Size: " << LSHSize(lsh) << endl;

    // 以後はLSHに格納されたデータを使うのでオリジナルはいらない
    cvReleaseMat(&objMat);

LSHの作成は、cvCreateMemoryLSH()という関数を使います。

    CvLSH* cvCreateMemoryLSH(int d, int n, int L, int k, int type, double r, int64 seed) 

dは、オリジナルデータの次元数です。今回はSURF特徴量なのでd=128ですね。nはハッシュテーブルのサイズです。これはかなり微妙なパラメータですが適当に100000にしてみました。値が大きいほどメモリは食いますが衝突がおきにくくなるので検索が高速になるようです。LとkはLSHの中核をなすパラメータです。Lはハッシュテーブルの数でkはハッシュキーの次元数です（k-NNのkとややこしいですが無関係）。d次元データは、k次元のハッシュキーに変換され、L個のハッシュテーブルに格納されます（実装ではL個のハッシュテーブルを作らず、サイズnのハッシュテーブルを1つ用意してデータの格納をL回繰り返しているようです）。LSHの詳しい解説は参考文献を見てください。今回の例では、L=5、k=64としています。何か最適のLとkの式が論文に書いてあるんですがその通りやったらうまくいきませんでした。今回はいろんなパラメータで試行錯誤Lとkを決めてます。ここら辺はもう少し追求したい。

LSHへデータを格納するのがcvLSHAdd()という関数です。

    void cvLSHAdd(CvLSH* lsh, const CvMat* data, CvMat* indices)

lshは、cvCreateMemoryLSHが返したLSHです。dataは、格納すべきデータです。1行に1データの行列（CvMat）に展開してからまとめてAddできます。今回の例だとdataは、局所特徴点の数xSURFの次元数（=128）の行列になります。kd-treeと違って、objMatのデータをコピーしてLSHに格納しているのでlshを作成した後はobjMatはいらないので解放してます。objMatが巨大だとobjMatとlshの両方で2倍のメモリを使うためデータ全部をobjMatに展開せず、少しずつディスクからロードしてcvLSHAdd()を繰り返したほうがメモリ効率はよさそうです。

        // クエリのキーポイントの特徴ベクトルをCvMatに展開
        CvMat* queryMat = cvCreateMat(queryDescriptors->total, DIM, CV_32FC1);
        CvSeqReader reader;
        float* ptr = queryMat->data.fl;
        cvStartReadSeq(queryDescriptors, &reader);
        for (int i = 0; i < queryDescriptors->total; i++) {
            float* descriptor = (float*)reader.ptr;
            CV_NEXT_SEQ_ELEM(reader.seq->elem_size, reader);
            memcpy(ptr, descriptor, DIM * sizeof(float));  // DIM次元の特徴ベクトルをコピー
            ptr += DIM;
        }

        // kd-treeで1-NNのキーポイントインデックスを検索
        int k = 1;  // k-NNのk
        CvMat* indices = cvCreateMat(queryKeypoints->total, k, CV_32SC1);   // 1-NNのインデックス
        CvMat* dists = cvCreateMat(queryKeypoints->total, k, CV_64FC1);     // その距離
        cvLSHQuery(lsh, queryMat, indices, dists, k, 100);

クエリの特徴ベクトル群を行列（CvMat）に展開して各特徴ベクトルの最近傍をkd-treeから検索する部分です。cvLSHQuery()という関数を使います。インタフェースは、kd-treeのcvFindFeatures()とほとんど同じです。LSHの作者があわせてくれたんですね。LSHのソースコードを読んでみましたが、ユークリッド距離空間のLSHを使っていました。実装は、Andoni氏のE2LSHと同じようです。

実行時間の比較

ためしに実験してみます。環境ですがCPUはCorei7 2.67GB、メモリ3GBです。前回の線形探索と同じCaltech101のサブセットを使います。

■kd-treeの結果
物体ID->物体名のハッシュを作成します ... OK
物体モデルデータベースをロードします ... OK
物体モデルデータベースをインデキシングします ... OK
物体モデルデータベースの物体数: 990
データベース中のキーポイント数: 321935
Loading Models Time = 78036.5ms
query? > dolphin_image_0001.jpg
../dataset/caltech101_10/dolphin_image_0001.jpg
クエリのキーポイント数: 183
識別結果: dolphin_image_0001.jpg
Recognition Time = 269.143ms
query? > butterfly_image_0001.jpg
../dataset/caltech101_10/butterfly_image_0001.jpg
クエリのキーポイント数: 423
識別結果: butterfly_image_0001.jpg
Recognition Time = 216.022ms
query? > 

■LSHの結果
物体ID->物体名のハッシュを作成します ... OK
物体モデルデータベースをロードします ... OK
物体モデルデータベースをインデキシングします ... OK
LSH Size: 321935
物体モデルデータベースの物体数: 990
データベース中のキーポイント数: 321935
Loading Models Time = 72197.6ms
query? > dolphin_image_0001.jpg
../dataset/caltech101_10/dolphin_image_0001.jpg
クエリのキーポイント数: 203
識別結果: dolphin_image_0001.jpg
Recognition Time = 51.9376ms
query? > butterfly_image_0001.jpg
../dataset/caltech101_10/butterfly_image_0001.jpg
クエリのキーポイント数: 453
識別結果: butterfly_image_0001.jpg
Recognition Time = 92.3854ms
query? > 

どちらもクエリと同じ画像を識別結果として返しているので正解ですね。線形探索と同様クエリ画像を少しくらい変形しても正しく認識できます。計算時間は、線形探索は23秒と52秒だったのでkd-treeとLSHは桁違いに速いことがわかります。kd-treeとLSHを比べるとLSHのほうが数倍早いです。もっと大きなデータだとさらに差がつきそうです。どちらも200MBくらいメモリを消費します（LSHはobjMatを解放するまでx2倍かかってます）。リアルタイム認識も普通のPCでもできました。

おわりに

ついこないだGoogleがGoogle Gogglesという物体認識を用いた検索システムを公開してました。これは予想ですが、今回作ったような高速特定物体認識が基盤になっているのではないかと思います。物体認識精度の改良や大規模データへの対応など実用化にあたっては困難な課題もありますが、簡単な物体認識システムなら素人でも3日で作れますね！

ちなみに東のエデンの実現はかなり大変だと思います。今回のシステムではマメシバや滝沢朗のリアルタイム認識は到底できません（笑）

参考文献

• lsh p-stable - LSHのわかりやすい解説資料、日本語では一番わかりやすかった
• LSH その1 - LSHの種類を解説した記事（原著論文へのリンクもあります）、上の解説資料を書いた方のブログ
• LSH Algorithm and Implementation (E2LSH) - Andoni氏による広く使われているユークリッド距離空間のLSH実装、LSHはハッシュキーがk次元ベクトルになりますがそれを普通のハッシュテーブルでどう実装すればよいかマニュアルに書かれてます
• Datar et al., Locality-Sensitive Hashing Scheme based on p-Stable Distributions（PDF）, Proceedings of the 20th Annual Symposium on Computational Geometry (2004) - 原著論文のPPT資料

このシリーズの記事一覧

• 3日で作る高速特定物体認識システム (1) 物体認識とは（2009/10/18）
• 3日で作る高速特定物体認識システム (2) SIFT特徴量の抽出（2009/10/24）
• 3日で作る高速特定物体認識システム (3) SURFの抽出（2009/10/30）
• 3日で作る高速特定物体認識システム (4) 特徴点のマッチング（2009/11/2）
• 3日で作る高速特定物体認識システム (5) 物体モデルデータベースの作成（2009/11/14）
• 3日で作る高速特定物体認識システム (6) 線形探索を用いた特定物体認識（2009/11/22）
• 3日で作る高速特定物体認識システム (7) 最近傍探索の高速化（2009/12/12）

■[OpenCV] 3日で作る高速特定物体認識システム (6) 線形探索を用いた特定物体認識 10:38

3日で作る高速特定物体認識システム (5) 物体モデルデータベースの作成（2009/11/14）のつづきです。

今回は、クエリを与えたときに物体モデルデータベースの中から一番似ている物体を探し出してクエリが何かを認識する部分（いわゆる、特定物体認識）を作ります。一番似ている物体を探し出す方法として、一番簡単な線形探索というアルゴリズムを使います。線形探索は遅いですが実装も簡単です。

キーポイントの最近傍点投票

クエリと一番似ている物体を検索するにはキーポイントの特徴量同士を比較します。各画像は数百から数千個の128次元の特徴ベクトルを持っているので単純に2つのベクトルを比較するという方法は使えません*1。代わりに最近傍点の投票方式を使います。クエリの各特徴ベクトルに一番近い（最近傍）特徴ベクトルを物体モデルデータベースから探し出しその特徴ベクトルを含む物体IDに対して投票します。たとえば、上の例では、クエリの特徴ベクトル(0)の最近傍を探して特徴ベクトル(1152)が見つかったとします。その特徴ベクトルを持つ物体はイルカなので物体ID=4に投票します（正の字で表してます）。クエリの残りの特徴ベクトルも同様です。例では物体ID=4のイルカへの投票が一番多いのでクエリの物体は物体ID=4のイルカと判定されます。

最近傍を探すのには線形探索を使います。これは、クエリの1つの特徴ベクトルの最近傍を求めるために物体モデルデータベースのすべての特徴ベクトルとの間でユークリッド距離を計算し、一番距離が小さいものを最近傍とします。なのでクエリ側に800個特徴ベクトルがあって、物体モデルデータベース側に300000個特徴ベクトルがあったら、投票が完了するまでに800x300000回距離計算が必要です。

クエリと物体モデルデータベースの画像が同じだったらまったく同じ特徴ベクトルが存在することになるので得票は1箇所に集中することになります。クエリ画像が少し変形していると得票は多少ばらつきますが一番得票数が多い物体を認識結果とすればOKです。特徴点のマッチング（2009/11/2）も参照してください。

プログラムリスト

ちょっと長いですが全プログラムリストです。

linear_recognition.cpp

#include <cv.h>
#include <highgui.h>
#include <iostream>
#include <fstream>
#include <map>

using namespace std;

const int DIM = 128;
const double DIST_THRESHOLD = 0.25;
const double VOTE_THRESHOLD = 50;
const int SURF_PARAM = 400;

const char* IMAGE_DIR = "caltech101_10";
const char* OBJID_FILE = "object.txt";
const char* DESC_FILE = "description.txt";

// プロトタイプ宣言
bool loadObjectId(const char *filename, map<int, string>& id2name);
bool loadDescription(const char *filename, vector<int> &labels, vector<int> &laplacians, CvMat* &objMat);
double euclideanDistance(float *vec, float *mvec, int length);
int searchNN(float *vec, int lap, vector<int> &labels, vector<int> &laplacians, CvMat* objMat);

int main(int argc, char** argv) {
    double tt = (double)cvGetTickCount();

    // 物体ID->物体ファイル名のハッシュを作成
    cout << "物体ID->物体名のハッシュを作成します ... " << flush;
    map<int, string> id2name;
    if (!loadObjectId(OBJID_FILE, id2name)) {
        cerr << "cannot load object id file" << endl;
        return 1;
    }
    cout << "OK" << endl;

    // キーポイントの特徴ベクトルをobjMat行列にロード
    cout << "物体モデルデータベースをロードします ... " << flush;
    vector<int> labels;     // キーポイントのラベル（objMatに対応）
    vector<int> laplacians;  // キーポイントのラプラシアン
    CvMat* objMat;           // 各行が物体のキーポイントの特徴ベクトル
    if (!loadDescription(DESC_FILE, labels, laplacians, objMat)) {
        cerr << "cannot load description file" << endl;
        return 1;
    }
    cout << "OK" << endl;

    cout << "物体モデルデータベースの物体数: " << id2name.size() << endl;
    cout << "データベース中のキーポイント数: " << objMat->rows << endl;
    tt = (double)cvGetTickCount() - tt;
    cout << "Loading Models Time = " << tt / (cvGetTickFrequency() * 1000.0) << "ms" << endl;

    while (1) {
        // クエリファイルの入力
        char input[1024];
        cout << "query? > ";
        cin >> input;

        char queryFile[1024];
        snprintf(queryFile, sizeof queryFile, "%s/%s", IMAGE_DIR, input);

        cout << queryFile << endl;

        tt = (double)cvGetTickCount();

        // クエリ画像をロード
        IplImage *queryImage = cvLoadImage(queryFile, CV_LOAD_IMAGE_GRAYSCALE);
        if (queryImage == NULL) {
            cerr << "cannot load image file: " << queryFile << endl;
            continue;
        }

        // クエリからSURF特徴量を抽出
        CvSeq *queryKeypoints = 0;
        CvSeq *queryDescriptors = 0;
        CvMemStorage *storage = cvCreateMemStorage(0);
        CvSURFParams params = cvSURFParams(SURF_PARAM, 1);
        cvExtractSURF(queryImage, 0, &queryKeypoints, &queryDescriptors, storage, params);
        cout << "クエリのキーポイント数: " << queryKeypoints->total << endl;

        // クエリの各キーポイントの1-NNの物体IDを検索して投票
        int numObjects = (int)id2name.size();  // データベース中の物体数
        int votes[numObjects];  // 各物体の集めた得票数
        for (int i = 0; i < numObjects; i++) {
            votes[i] = 0;
        }
        for (int i = 0; i < queryDescriptors->total; i++) {
            CvSURFPoint *p = (CvSURFPoint *)cvGetSeqElem(queryKeypoints, i);
            float *vec = (float *)cvGetSeqElem(queryDescriptors, i);
            int lap = p->laplacian;
            int nnId = searchNN(vec, lap, labels, laplacians, objMat);
            votes[nnId]++;
        }

        for (int i = 0; i < numObjects; i++) {
            cout << votes[i] << endl;
        }

        // 投票数が最大の物体IDを求める
        int maxId = -1;
        int maxVal = -1;
        for (int i = 0; i < numObjects; i++) {
            if (votes[i] > maxVal) {
                maxId = i;
                maxVal = votes[i];
            }
        }

        // 物体IDを物体ファイル名に変換
        string name = id2name[maxId];
        cout << "識別結果: " << name << endl;

        tt = (double)cvGetTickCount() - tt;
        cout << "Recognition Time = " << tt / (cvGetTickFrequency() * 1000.0) << "ms" << endl;

        // 後始末
        cvReleaseImage(&queryImage);
        cvClearSeq(queryKeypoints);
        cvClearSeq(queryDescriptors);
        cvReleaseMemStorage(&storage);
        cvDestroyAllWindows();
    }

    // 後始末
    cvReleaseMat(&objMat);

    return 0;
}

/**
 * 物体ID->物体名のmapを作成して返す
 *
 * @param[in]  filename  物体ID->物体名の対応を格納したファイル
 * @param[out] id2name   物体ID->物体名のmap
 *
 * @return 成功ならtrue、失敗ならfalse
 */
bool loadObjectId(const char *filename, map<int, string>& id2name) {
    // 物体IDと物体名を格納したファイルを開く
    ifstream objFile(filename);
    if (objFile.fail()) {
        cerr << "cannot open file: " << filename << endl;
        return false;
    }

    // 1行ずつ読み込み、物体ID->物体名のmapを作成
    string line;
    while (getline(objFile, line, '\n')) {
        // タブで分割した文字列をldataへ格納
        vector<string> ldata;
        istringstream ss(line);
        string s;
        while (getline(ss, s, '\t')) {
            ldata.push_back(s);
        }

        // 物体IDと物体名を抽出してmapへ格納
        int objId = atol(ldata[0].c_str());
        string objName = ldata[1];
        id2name.insert(map<int, string>::value_type(objId, objName));
    }

    // 後始末
    objFile.close();

    return true;
}

/**
 * キーポイントのラベル（抽出元の物体ID）とラプラシアンと特徴ベクトルをロードしlabelsとobjMatへ格納
 *
 * @param[in]  filename         特徴ベクトルを格納したファイル
 * @param[out] labels           特徴ベクトル抽出元の物体ID
 * @param[out] laplacianss      特徴ベクトルのラプラシアン
 * @param[out] objMat           特徴量を格納した行列（各行に1つの特徴ベクトル）
 *
 * @return 成功ならtrue、失敗ならfalse
 */
bool loadDescription(const char *filename, vector<int> &labels, vector<int> &laplacians, CvMat* &objMat) {
    // 物体IDと特徴ベクトルを格納したファイルを開く
    ifstream descFile(filename);
    if (descFile.fail()) {
        cerr << "cannot open file: " << filename << endl;
        return false;
    }

    // 行列のサイズを決定するためキーポイントの総数をカウント
    int numKeypoints = 0;
    string line;
    while (getline(descFile, line, '\n')) {
        numKeypoints++;
    }
    objMat = cvCreateMat(numKeypoints, DIM, CV_32FC1);

    // ファイルポインタを先頭に戻す
    descFile.clear();
    descFile.seekg(0);

    // データを読み込んで行列へ格納
    int cur = 0;
    while (getline(descFile, line, '\n')) {
        // タブで分割した文字列をldataへ格納
        vector<string> ldata;
        istringstream ss(line);
        string s;
        while (getline(ss, s, '\t')) {
            ldata.push_back(s);
        }
        // 物体IDを取り出して特徴ベクトルのラベルとする
        int objId = atol(ldata[0].c_str());
        labels.push_back(objId);
        // ラプラシアンを取り出して格納
        int laplacian = atoi(ldata[1].c_str());
        laplacians.push_back(laplacian);
        // DIM次元ベクトルの要素を行列へ格納
        for (int j = 0; j < DIM; j++) {
            float val = atof(ldata[j+2].c_str());  // 特徴ベクトルはldata[2]から
            CV_MAT_ELEM(*objMat, float, cur, j) = val;
        }
        cur++;
    }

    descFile.close();

    return true;
}

/**
 * 2つのベクトルのユークリッド距離を計算して返す
 *
 * @param[in] vec     ベクトル1の配列
 * @param[in] mvec    ベクトル2の配列
 * @param[in] length  ベクトルの長さ
 *
 * @return ユークリッド距離
 */
double euclideanDistance(float *vec, float *mvec, int length) {
    double dist = 0.0;

    for (int i = 0; i < length; i++) {
        dist += (vec[i] - mvec[i]) * (vec[i] - mvec[i]);
    }
    dist = sqrt(dist);

    return dist;
}

/**
 * クエリのキーポイントの1-NNキーポイントを物体モデルデータベースから探してその物体IDを返す
 *
 * @param[in] vec          クエリキーポイントの特徴ベクトル
 * @param[in] lap          クエリキーポイントのラプラシアン
 * @param[in] labels       物体モデルデータベースの各キーポイントの物体ID
 * @param[in] laplacians   物体モデルデータベースの各キーポイントのラプラシアン
 * @param[in] objMat       物体モデルデータベースの各キーポイントの特徴ベクトル
 *
 * @return 指定したキーポイントにもっとも近いキーポイントの物体ID
 */
int searchNN(float *vec, int lap, vector<int> &labels, vector<int> &laplacians, CvMat* objMat) {
    int neighborObj = -1;
    double minDist = 1e6;
    for (int i = 0; i < objMat->rows; i++) {
        // クエリのキーポイントとラプラシアンが異なるキーポイントは無視
        if (lap != laplacians[i]) {
            continue;
        }
        // i行目の先頭データのポインタ（ここから128個がi行目の特徴ベクトル）
        float *mvec = (float *)(objMat->data.fl + i * DIM);
        double d = euclideanDistance(vec, mvec, DIM);
        if (d < minDist) {
            neighborObj = labels[i];  // NNキーポイントの物体IDを更新
            minDist = d;           // 最小の距離を更新
        }
    }
    return neighborObj;
}

loadObjectId()とloadDescription()は、物体モデルデータベースの作成（2009/11/14）で作ったobject.txtとdescripion.txtをロードする関数です。両ファイルと画像フォルダはプログラムと同じフォルダに置いてください。loadObjectId()は、物体ID->ファイル名のmapを作成します。認識結果を表示するときに「クエリの物体IDは3です。」とか言われても分かりにくいので名前を表示するために使います。loadDescription()は、特徴量ファイルを読み込んで特徴ベクトルの物体IDをlabelsにラプラシアンをlaplaciansに特徴ベクトルをobjMatにロードします。searchNN()は、vecとlapにクエリの特徴ベクトルとラプラシアンを指定すると物体モデルデータベースを線形探索して一番近い特徴ベクトルを検索し、その物体IDを返します。main()ではクエリ画像のファイル名を与えるとその画像から特徴量を抽出、クエリの各特徴ベクトルの最近傍を検索して物体IDの得票数を計算し、得票数がもっとも多い物体IDのファイル名を認識結果として表示してます。

実行例

以下、プログラムの実行例です。クエリに

• dolphin_image_0001.jpg
• butterfly_image_0001.jpg

を与えてみました。170MBくらいメモリを使うので少ない人は気をつけてください。もっとサイズが小さいデータセットで試してみるとよいかも。

物体ID->物体名のハッシュを作成します ... OK
物体モデルデータベースをロードします ... OK
物体モデルデータベースの物体数: 990
データベース中のキーポイント数: 321935
Loading Models Time = 64689.6ms
query? > dolphin_image_0001.jpg
../dataset/caltech101_10/dolphin_image_0001.jpg
クエリのキーポイント数: 203
識別結果: dolphin_image_0001.jpg
Recognition Time = 23301.2ms
query? > butterfly_image_0001.jpg
../dataset/caltech101_10/butterfly_image_0001.jpg
クエリのキーポイント数: 453
識別結果: butterfly_image_0001.jpg
Recognition Time = 52065.9ms
query? > 

ちゃんと認識結果がクエリと同じファイルになることがわかると思います。でもこれは物体モデルデータベースに登録された画像と同じ画像をクエリにしているので当たり前の結果です。実はクエリ画像を少し変形して（回転、縮小、拡大、隠蔽など）与えてもちゃんと認識できます。線形探索なので計算はかなり遅いです。クエリを与えてから認識結果を返すまで23秒と52秒かかってます。こんなに遅くては到底リアルタイムで認識させられませんねー。

これでやっと3日で作る高速特定物体認識システムのうち1日目の課題が終わたってことですね。今回使った線形探索は遅すぎるので、物体モデルデータベースの局所特徴量をkd-treeやLocality Sensitive Hashing（LSH）でインデキシングして最近傍探索を高速化することを考えていきます。線形探索とは比べ物にならないくらい計算が速くなります。

3日で作る高速特定物体認識システム (7) 最近傍探索の高速化（2009/12/12）へつづきます。次回でこのレポート最終回の予定です。

このシリーズの記事一覧

• 3日で作る高速特定物体認識システム (1) 物体認識とは（2009/10/18）
• 3日で作る高速特定物体認識システム (2) SIFT特徴量の抽出（2009/10/24）
• 3日で作る高速特定物体認識システム (3) SURFの抽出（2009/10/30）
• 3日で作る高速特定物体認識システム (4) 特徴点のマッチング（2009/11/2）
• 3日で作る高速特定物体認識システム (5) 物体モデルデータベースの作成（2009/11/14）
• 3日で作る高速特定物体認識システム (6) 線形探索を用いた特定物体認識（2009/11/22）
• 3日で作る高速特定物体認識システム (7) 最近傍探索の高速化（2009/12/12）

■[OpenCV] 3日で作る高速特定物体認識システム (5) 物体モデルデータベースの作成 13:17

3日で作る高速特定物体認識システム (4) 特徴点のマッチング（2009/11/2）のつづきです。

今回は、高速特定物体認識システムの物体モデルデータベースを作成してみます。クエリとして与えた画像が何かはこの物体モデルデータベースを検索することで認識できるようになります。特定物体認識では、クエリとして与えた画像が物体モデルデータベースにあらかじめ登録されていないとまず認識できません（局所特徴量を用いるので多少のスケール変化や回転などには対応できますが）。忘れてしまった方は物体認識システムの構成（2009/10/18）を参照してください。

画像データセットの準備

自分で撮影した画像を登録できるようにしてもよかったのですが、大量の画像を撮影するのは面倒なのでおなじみのCaltech101の画像セットを用いて物体モデルを作成します。

類似画像検索システムを作ろう（2009/10/3）の画像データセットを用意を参考にcaltech101フォルダに画像ファイルだけコピーしたデータを用意します。少し多すぎるので各カテゴリから10枚ずつ画像を抜き出して（IDが0001から0010まで）caltech101_10というデータを作ります。caltech101_10というフォルダを作成して以下のスクリプトを実行すると990枚の画像データができます。

#coding:utf-8
import os
import shutil

IMAGE_DIR = "caltech101"
OUT_DIR = "caltech101_10"

for file in os.listdir(IMAGE_DIR):
    id = int(file[-8:-4])
    if id <= 10:
        print file
        shutil.copy("%s/%s" % (IMAGE_DIR, file), "%s/%s" % (OUT_DIR, file))

物体モデルデータベースの作成

caltech101_10という990枚の画像から特徴量を抽出して物体モデルデータベースに格納していきます。各画像から抽出したキーポイントの特徴量は物体IDとひもづけて保存する必要があります。そこで、物体IDと特徴量を保存するテーブルが必要になります。

物体IDテーブル：物体ID 画像ファイル名
特徴量テーブル：物体ID ラプラシアン キーポイントの特徴量

SQLiteなどでデータベースを作ってもよいのですが、コードが複雑になるのでここでは簡単のためファイルにテキスト形式でデータを保存してみました。物体IDファイル（object.txt）と特徴量ファイル（description.txt）です。

物体IDファイルは以下のようなタブ区切りの書式です。登録した順番に物体IDが自動で割り当てられます。

0    accordion_image_0001.jpg
1    ant_image_0001.jpg
2    bonsai_image_0001.jpg
3    butterfly_image_0001.jpg
4    dolphin_image_0001.jpg

特徴量ファイルは以下のようなタブ区切りの書式です。すべての画像の特徴量は1つのファイルに保存しています。各物体は複数のキーポイントを持つ点に注意が必要です。なのでデータベースを使う場合は物体IDを主キーにできません。

0    1    0.03    0.05    ....   # 物体ID=0の1つ目のキーポイント
0    1    0.28    0.14    ....   # 物体ID=0の2つ目のキーポイント
0   -1    0.31    0.01    ....   # 物体ID=0の3つ目のキーポイント
....
1   -1    0.24    0.31    ....   # 物体ID=1の1つ目のキーポイント
1    1    0.15    0.22    ....   # 物体ID=1の2つ目のキーポイント
1    1    0.01    0.03    ....   # 物体ID=1の3つ目のキーポイント
....

caltech101_10の全画像ファイルから上記の物体IDファイルと特徴量ファイルを抽出するプログラムです。局所特徴量にはSURFを使っています。詳しくは、SURFの抽出（2009/10/30）を参照してください。ディレクトリのファイルスキャンにopendir()を使っています。これは、MinGWのg++では使えますが、もしかしたらVisual C++では使えないかも・・・もしエラーが出たらここらを参考に書き換えてください。

#include <cv.h>
#include <highgui.h>
#include <iostream>
#include <fstream>
#include <sys/types.h>
#include <dirent.h>

using namespace std;

const char *IMAGE_DIR = "caltech101_10";
const char *OBJ_FILE = "object.txt";        // 物体ID格納ファイル
const char *DESC_FILE = "description.txt";  // 特徴量格納ファイル

const double SURF_PARAM = 400;  // SURFのパラメータ
const int DIM = 128;            // SURF特徴量の次元数

/**
 * SURF特徴量を抽出する
 *
 * @param[in]  filename         画像ファイル名
 * @param[out] imageKeypoints   キーポイント（出力のため参照渡し）
 * @param[out] imageDescriptors 各キーポイントのSURF特徴量（出力のため参照渡し）
 * @param[out] storage          Memory Storage（出力のため参照渡し）
 *
 * @return 成功なら0、失敗なら1
 */
int extractSURF(char *filename, CvSeq* &imageKeypoints, CvSeq* &imageDescriptors, CvMemStorage* &storage) {
    // グレースケールで画像をロードする
    IplImage *img = cvLoadImage(filename, CV_LOAD_IMAGE_GRAYSCALE);
    if (img == NULL) {
        cerr << "cannot load image file: " << filename << endl;
        return 1;
    }

    storage = cvCreateMemStorage(0);
    CvSURFParams params = cvSURFParams(SURF_PARAM, 1);
    cvExtractSURF(img, 0, &imageKeypoints, &imageDescriptors, storage, params);

    return 0;
}

/**
 * 物体モデルをファイルに保存する
 *
 * @param[in]   objId              オブジェクトID
 * @param[in]   filename           画像ファイル名
 * @param[in]   imageKeypoints     キーポイント
 * @param[in]   imageDescriptors   各キーポイントの特徴量
 * @param[in]   objFile            物体IDファイルのハンドラ
 * @param[in]   descFile           特徴量ファイルのハンドラ
 *
 * @return 成功なら0、失敗なら1
 */
int saveFile(int objId, char *filename, CvSeq* imageKeypoints, CvSeq* imageDescriptors, ofstream& objFile, ofstream& descFile) {
    cout << objId << " " << filename << " " << imageDescriptors->total << endl;

    // 物体IDファイルへ登録
    objFile << objId << "\t" << filename << endl;

    // オブジェクトID, ラプラシアン, 128個の数字をタブ区切りで出力
    for (int i = 0; i < imageDescriptors->total; i++) {  // 各キーポイントの特徴量に対し
        // オブジェクトID
        descFile << objId << "\t";

        // 特徴点のラプラシアン（SURF特徴量ではベクトルの比較時に使用）
        const CvSURFPoint* kp = (const CvSURFPoint*)cvGetSeqElem(imageKeypoints, i);
        int laplacian = kp->laplacian;
        descFile << laplacian << "\t";

        // 128次元ベクトル
        const float *descriptor = (const float *)cvGetSeqElem(imageDescriptors, i);
        for (int d = 0; d < DIM; d++) {
            descFile << descriptor[d] << "\t";
        }

        descFile << endl;
    }

    return 0;
}

int main(int argc, char **argv) {
    int ret;

    // 物体IDファイルを開く
    ofstream objFile(OBJ_FILE);
    if (objFile.fail()) {
        cerr << "cannot open file: " << OBJ_FILE << endl;
        return 1;
    }

    // 特徴量ファイルを開く
    ofstream descFile(DESC_FILE);
    if (descFile.fail()) {
        cerr << "cannot open file: " << DESC_FILE << endl;
        return 1;
    }

    // IMAGE_DIRの画像ファイル名を走査
    DIR *dp = opendir(IMAGE_DIR);
    if (dp == NULL) {
        cerr << "cannot open directory: " << IMAGE_DIR << endl;
        return 1;
    }

    int objId = 0;  // オブジェクトID
    struct dirent *entry;
    while (1) {
        entry = readdir(dp);

        if (entry == NULL) {
            break;
        }

        // .と..は無視する
        if (strncmp(entry->d_name, ".", 1) == 0 || strncmp(entry->d_name, "..", 2) == 0) {
            continue;
        }

        char *filename = entry->d_name;

        // SURFを抽出
        char buf[1024];
        snprintf(buf, sizeof buf, "%s/%s", IMAGE_DIR, filename);
        CvSeq *imageKeypoints = 0;
        CvSeq *imageDescriptors = 0;
        CvMemStorage *storage = 0;
        ret = extractSURF(buf, imageKeypoints, imageDescriptors, storage);
        if (ret != 0) {
            cerr << "cannot extract surf description" << endl;
            return 1;
        }

        // ファイルへ出力
        ret = saveFile(objId, filename, imageKeypoints, imageDescriptors, objFile, descFile);
        if (ret != 0) {
            cerr << "cannot save surf description" << endl;
            return 1;
        }

        // 後始末
        cvClearSeq(imageKeypoints);
        cvClearSeq(imageDescriptors);
        cvReleaseMemStorage(&storage);

        objId++;
    }

    objFile.close();
    descFile.close();
    closedir(dp);

    return 0;
}

上記のプログラムで990枚の画像から物体IDファイルと特徴量ファイルの作成にはうちの環境で5分くらいかかりました。description.txtはかなりでかくて420MBもあります。うかつにメモ帳とかで開かないようにしないと（笑）サイズがでかいのは浮動小数点数をテキストで保存してるからですが、バイナリで格納すればもっと小さくなりますね。大規模な画像を対象にする場合は工夫が必要です。990枚の画像の合計キーポイント数は321932個ありました。1枚あたり平均325個という計算です。Caltech101の画像はサイズが小さいので少ないです。

物体モデルデータベースを特定物体認識にどう使うか？

今回作成したobject.txtとdescription.txtをここでは物体モデルデータベースと呼びます。これを特定物体認識システムにどう使うかですが下の図のような感じです。

今回は990枚の画像を格納した図右側の物体モデルデータベースを作成しました。特定物体認識ではクエリ画像（左のイルカ）を与えてこれが物体モデルデータベースのどの画像かを言い当てるという処理です。上の例の場合、オブジェクトID=4と認識できれば大正解です。

この認識方法は、すでにご紹介している特徴点のマッチング（2009/11/2）とほとんど同じ方法です。特徴点のマッチングでは画像対画像でしたが、特定物体認識では画像対物体モデルデータベースで比較するという違いがあるだけです。

物体モデルデータベースには大量のキーポイントが含まれているため特徴点のマッチングより高速な最近傍計算手法が必要になります。特徴点のマッチングでは画像1のキーポイント数が800、画像2のキーポイント数が800だとすると全部で800x800=640000回程度の距離計算で画像1の各キーポイントの最近傍点が見つかります。ここで画像2を物体モデルデータベースに置き換えると全部で800x321932=257545600回もの距離計算になります。効率の悪い線形探索では時間がかかりすぎるわけです。とはいうもののどれくらいかかるのか確かめるため次回は線形探索を用いた特定物体認識を行ってみます。

3日で作る高速特定物体認識システム (6) 線形探索を用いた特定物体認識へつづきます。

このシリーズの記事一覧

• 3日で作る高速特定物体認識システム (1) 物体認識とは（2009/10/18）
• 3日で作る高速特定物体認識システム (2) SIFT特徴量の抽出（2009/10/24）
• 3日で作る高速特定物体認識システム (3) SURFの抽出（2009/10/30）
• 3日で作る高速特定物体認識システム (4) 特徴点のマッチング（2009/11/2）
• 3日で作る高速特定物体認識システム (5) 物体モデルデータベースの作成（2009/11/14）
• 3日で作る高速特定物体認識システム (6) 線形探索を用いた特定物体認識（2009/11/22）
• 3日で作る高速特定物体認識システム (7) 最近傍探索の高速化（2009/12/12）

■[OpenCV] 3日で作る高速特定物体認識システム (4) 特徴点のマッチング 22:09

3日で作る高速特定物体認識システム (3) SURFの抽出（2009/10/30）のつづき。

画像からSIFTや SURFといった局所特徴量を抽出できるようになったのでここらでそれを応用してみます。特徴点のマッチングを取ることで2つの画像間で対応する場所を求められるようになります。下の例のような感じです。下の図で2つのキーポイント間にひいた直線は、両端のキーポイントの特徴ベクトルが似ている（距離が小さい）ことを表しています。

以下、プログラムです。

keypoint_matching.exe [画像1のファイル名] [画像2のファイル名]

のように2つの画像ファイルを入力として与えると上のようにマッチング画像が表示されます。

#include <cv.h>
#include <highgui.h>
#include <iostream>
#include <vector>
#include <cmath>

using namespace std;

const int DIM_VECTOR = 128;    // 128次元ベクトル
const double THRESHOLD = 0.3;  // 線をつなぐ距離の閾値

/**
 * 2つのベクトルのユークリッド距離を計算して返す
 *
 * @param[in] vec     ベクトル1の配列
 * @param[in] mvec    ベクトル2の配列
 * @param[in] length  ベクトルの長さ
 *
 * @return ユークリッド距離
 */
double euclidDistance(float* vec1, float* vec2, int length) {
    double sum = 0.0;
    for (int i = 0; i < length; i++) {
        sum += (vec1[i] - vec2[i]) * (vec1[i] - vec2[i]);
    }
    return sqrt(sum);
}

/**
 * 最近傍点を検索
 *
 * @param[in]   vec             特徴ベクトル
 * @param[in]   laplacian       ラプラシアン
 * @param[in]   keypoints       キーポイントの集合（この中から最近傍点を検索）
 * @param[in]   descriptors     特徴ベクトルの集合
 *
 * @return 最近傍点のインデックス（見つからないときは-1）
 */
int nearestNeighbor(float* vec, int laplacian, CvSeq* keypoints, CvSeq* descriptors) {
    int neighbor = -1;
    double minDist = 1e6;

    for (int i = 0; i < descriptors->total; i++) {
        CvSURFPoint* pt = (CvSURFPoint*)cvGetSeqElem(keypoints, i);
        // ラプラシアンが異なるキーポイントは無視
        if (laplacian != pt->laplacian) continue;
        float* v = (float*)cvGetSeqElem(descriptors, i);
        double d = euclidDistance(vec, v, DIM_VECTOR);
        // より近い点があったら置き換える
        if (d < minDist) {
            minDist = d;
            neighbor = i;
        }
    }

    // 最近傍点でも距離が閾値以上だったら無視する
    if (minDist < THRESHOLD) {
        return neighbor;
    }

    // 最近傍点が見つからない場合
    return -1;
}

/**
 * 画像1のキーポイントと近い画像2のキーポイントを組にして返す
 *
 * @param[in]  keypoints1       画像1のキーポイント
 * @param[in]  descriptors1     画像1の特徴ベクトル
 * @param[in]  keypoints2       画像2のキーポイント
 * @param[in]  descriptors2     画像2の特徴ベクトル
 * @param[out] ptpairs          類似キーポイントインデックスの列（2つ1組）
 *
 * @return なし
 */
void findPairs(CvSeq* keypoints1, CvSeq* descriptors1,
                CvSeq* keypoints2, CvSeq* descriptors2,
                vector<int>& ptpairs) {
    ptpairs.clear();
    // 画像1の各キーポイントに関して最近傍点を検索
    for (int i = 0; i < descriptors1->total; i++) {
        CvSURFPoint* pt1 = (CvSURFPoint*)cvGetSeqElem(keypoints1, i);
        float* desc1 = (float*)cvGetSeqElem(descriptors1, i);
        // 最近傍点を検索
        int nn = nearestNeighbor(desc1, pt1->laplacian, keypoints2, descriptors2);
        // 見つかったら画像1のインデックスと画像2のインデックスを順番に登録
        if (nn >= 0) {
            ptpairs.push_back(i);
            ptpairs.push_back(nn);
        }
    }
}

int main(int argc, char** argv) {
    const char* filename1 = argc == 3 ? argv[1] : "image/dolphin_image_0001.jpg";
    const char* filename2 = argc == 3 ? argv[2] : "image/dolphin_image_0001.jpg";

    cvNamedWindow("Keypoint Matching");

    // 画像はグレースケールでロード
    IplImage* grayImage1 = cvLoadImage(filename1, CV_LOAD_IMAGE_GRAYSCALE);
    IplImage* grayImage2 = cvLoadImage(filename2, CV_LOAD_IMAGE_GRAYSCALE);
    if (!grayImage1 || !grayImage2) {
        cerr << "cannot find image file" << endl;
        exit(-1);
    }

    // 描画用にカラーでロード
    IplImage* colorImage1 = cvLoadImage(filename1, CV_LOAD_IMAGE_ANYCOLOR|CV_LOAD_IMAGE_ANYDEPTH);
    IplImage* colorImage2 = cvLoadImage(filename2, CV_LOAD_IMAGE_ANYCOLOR|CV_LOAD_IMAGE_ANYDEPTH);

    // マッチング用の画像を描画
    CvSize sz = cvSize(colorImage1->width + colorImage2->width, colorImage1->height + colorImage2->height);
    IplImage* matchingImage = cvCreateImage(sz, IPL_DEPTH_8U, 3);

    // 画像1を描画
    cvSetImageROI(matchingImage, cvRect(0, 0, colorImage1->width, colorImage1->height));
    cvCopy(colorImage1, matchingImage);
    // 画像2を描画
    cvSetImageROI(matchingImage, cvRect(colorImage1->width, colorImage1->height, colorImage2->width, colorImage2->height));
    cvCopy(colorImage2, matchingImage);
    cvResetImageROI(matchingImage);

    CvSeq *keypoints1 = 0, *descriptors1 = 0;
    CvSeq *keypoints2 = 0, *descriptors2 = 0;
    CvMemStorage* storage = cvCreateMemStorage(0);
    CvSURFParams params = cvSURFParams(500, 1);

    // SURFを抽出
    cvExtractSURF(grayImage1, 0, &keypoints1, &descriptors1, storage, params);
    cvExtractSURF(grayImage2, 0, &keypoints2, &descriptors2, storage, params);

    // 特徴ベクトルの類似度が高いキーポイント同士を線でつなぐ
    vector<int> ptpairs;  // keypointsのインデックスが2つずつ組になるように格納
    findPairs(keypoints1, descriptors1, keypoints2, descriptors2, ptpairs);

    // 2つずつ組にして線を引く
    for (int i = 0; i < (int)ptpairs.size(); i += 2) {
        CvSURFPoint* pt1 = (CvSURFPoint*)cvGetSeqElem(keypoints1, ptpairs[i]);     // 画像1のキーポイント
        CvSURFPoint* pt2 = (CvSURFPoint*)cvGetSeqElem(keypoints2, ptpairs[i + 1]); // 画像2のキーポイント
        CvPoint from = cvPointFrom32f(pt1->pt);
        CvPoint to = cvPoint(cvRound(colorImage1->width + pt2->pt.x), cvRound(colorImage1->height + pt2->pt.y));
        cvLine(matchingImage, from, to, cvScalar(0, 255, 255));
    }

    // キーポイントマッチングした画像を描画
    cvShowImage("Keypoint Matching", matchingImage);
    cvWaitKey(0);

    // 後始末
    cvReleaseImage(&grayImage1);
    cvReleaseImage(&grayImage2);
    cvReleaseImage(&colorImage1);
    cvReleaseImage(&colorImage2);
    cvClearSeq(keypoints1);
    cvClearSeq(descriptors1);
    cvClearSeq(keypoints2);
    cvClearSeq(descriptors2);
    cvReleaseMemStorage(&storage);

    return 0;
}

画像1（左上の画像）の各キーポイントについて画像2（右下の画像）から最近傍点を検索しています。検索の方法は線形探索で、すべてのキーポイント間で距離を計算して一番近いものを選んでます。そのため、画像1のキーポイント数が800、画像2のキーポイント数が800だとすると全部で800x800=640000回もの距離計算が必要になります。線形探索は非常に効率が悪いので後でさらに高速な検索方法を導入予定です。ここで、SURFはラプラシアンが一致しない点同士は似ていないと判定してよいそうです。なぜかは勉強中です・・・

キーポイント間の距離は特徴ベクトル（128次元）のユークリッド距離（L2ノルム）で測っています。そして、ユークリッド距離がTHRESHOLD（=0.3）以下の場合、対応する点と判断して直線を引いています。実際は、距離の比較なのでユークリッド距離の2乗を使っても問題ありません（閾値は変える必要あり）。大量に距離計算が必要なのでsqrt()を使わずint型で計算したほうが高速かも。

局所特徴量の頑健性の検証

SIFTやSURFは画像のスケール変化、平行移動、回転、隠蔽（オクルージョン）に対して頑健と言われてます。そんなわけで画像1（左上の画像）を変形してキーポイントのマッチングができるか試してみました。画像をクリックで拡大できます。

画像のサイズ変更、回転、隠蔽は下のようなコードでできます。隠蔽は左半分を白く塗りつぶしています。詳しくは下のリンクを参照してください。

• 画像のサイズ変更
• 画像の回転
• 図形の描画

    // 縮小
    IplImage* dst = cvCreateImage(cvSize(src->width/2, src->height/2), src->depth, src->nChannels);
    cvResize(src, dst);

    // 隠蔽（左半分を塗りつぶす）
    IplImage* dst = cvCloneImage(src);
    cvRectangle(dst, cvPoint(0, 0), cvPoint(dst->width/2, dst->height), cvScalar(255,255,255), CV_FILLED);

    // 回転
    IplImage* dst = cvCloneImage(src);
    int angle = 45;  // 回転角度
    float m[6];
    CvMat M;
    m[0] = (float)(cos(angle * CV_PI / 180));
    m[1] = (float)(-sin(angle * CV_PI / 180));
    m[2] = src->width * 0.5;
    m[3] = -m[1];
    m[4] = m[0];
    m[5] = src->height * 0.5;
    cvInitMatHeader(&M, 2, 3, CV_32FC1, m, CV_AUTOSTEP);
    cvGetQuadrangleSubPix(src, dst, &M);

特定物体認識への応用

キーポイントのマッチングは物体画像があるシーンの一部にある場合でも可能です。たとえば、画像2（右下の画像）をいろいろな物体を結合した画像にしてみました。線を引くユークリッド距離の閾値を0.2にしています。ちゃんとイルカがどこにあるか認識してますね。現状では、画像2（右下の画像）はたくさんの物体画像をまとめた一枚画ですが、後ほど物体モデルデータベース、大量の物体画像のキーポイントと特徴ベクトルを格納したデータベースに置き換える予定です。

次の例は、ユークリッド距離の閾値を0.3にあげて少し条件を緩くしてみました。こうすると線はたくさん引かれますがその分間違えも大きくなってます。

別の例です。これもユークリッド距離の閾値は0.3です。間違いは多いですが、画像1と同じピアノの画像に線が集中しているのはわかります。実際に特徴点のマッチングを特定物体認識に使う場合は間違えが多少あっても許されます。下のように正解の画像に線が集中していれば（得票数が多ければ）きちんと認識できるわけです。この多少間違えがあっても認識できるという性質を利用して認識を高速化させる方法も提案されています。

3日で作る高速特定物体認識システム (5) 物体モデルデータベースの作成（2009/11/14）へ続きます。

このシリーズの記事一覧

• 3日で作る高速特定物体認識システム (1) 物体認識とは（2009/10/18）
• 3日で作る高速特定物体認識システム (2) SIFT特徴量の抽出（2009/10/24）
• 3日で作る高速特定物体認識システム (3) SURFの抽出（2009/10/30）
• 3日で作る高速特定物体認識システム (4) 特徴点のマッチング（2009/11/2）
• 3日で作る高速特定物体認識システム (5) 物体モデルデータベースの作成（2009/11/14）
• 3日で作る高速特定物体認識システム (6) 線形探索を用いた特定物体認識（2009/11/22）
• 3日で作る高速特定物体認識システム (7) 最近傍探索の高速化（2009/12/12）

■[OpenCV] 3日で作る高速特定物体認識システム (3) SURFの抽出 21:20

3日で作る高速特定物体認識システム (2) SIFT特徴量の抽出（2009/10/24）の続きです。あっ、3日経っちゃいました。

今回は、SIFTとは別の局所特徴量であるSURF（Speeded Up Robust Features）を抽出してみます。SURFのFはFeaturesなのでSURF特徴量とは言わないのかな？SIFTとは抽出方法は違いますが、画像からキーポイントと特徴ベクトルを抽出する点では同じです。抽出速度はSIFTより数倍高速だそうですが、精度は多少落ちるとのこと。リアルタイム処理したいときはこっちのほうがよさそうです。また、OpenCVにもすでに実装されています。SURFの詳しいアルゴリズムは後で論文を読むとしてとりあえず試してみます。

画像からSURFを抽出する

以下のプログラムは、画像からSURFを抽出して特徴点を描画し、特徴量をファイルへ格納するプログラムです。このプログラムはOpenCVについているfind_obj.cpp（C:\OpenCV2.0\samples\cにあります）を参考にしています。このプログラムからSURFを抽出する部分だけ取り出してみました。

surf.exe [画像ファイル名] [SURFファイル名]

という感じで使います。

#include <cv.h>
#include <highgui.h>
#include <iostream>
#include <fstream>

using namespace std;

const int DIM_VECTOR = 128;  // 128次元ベクトル

/**
 * SURF情報をファイルに出力
 * @param[in]   filename            SURFを保存するファイル名
 * @param[in]   imageKeypoints      SURFキーポイント情報
 * @param[in]   imageDescriptors    SURF特徴ベクトル情報
 * @return なし
 */
void writeSURF(const char* filename, CvSeq* imageKeypoints, CvSeq* imageDescriptors) {
    fstream fout;

    fout.open(filename, ios::out);
    if (!fout.is_open()) {
        cerr << "cannot open file: " << filename << endl;
        return;
    }

    // 1行目はキーポイント数と特徴量の次元数を書き込む
    fout << imageKeypoints->total << ' ' << DIM_VECTOR << endl;

    // 2行目からキーポイント情報と特徴ベクトルを書き込む
    for (int i = 0; i < imageKeypoints->total; i++) {
        CvSURFPoint* point = (CvSURFPoint*)cvGetSeqElem(imageKeypoints, i);
        float* descriptor = (float*)cvGetSeqElem(imageDescriptors, i);
        // キーポイント情報（X座標, Y座標, サイズ, ラプラシアン）を書き込む
        fout << point->pt.x << ' ' << point->pt.y << ' ' << point->size << ' ' << point->laplacian << ' ';
        // 特徴ベクトルを書き込む
        for (int j = 0; j < DIM_VECTOR; j++) {
            fout << descriptor[j] << ' ';
        }
        fout << endl;
    }

    fout.close();
}

int main(int argc, char** argv) {
    const char* imageFile = argc == 3 ? argv[1] : "image/accordion_image_0001.jpg";
    const char* surfFile  = argc == 3 ? argv[2] : "image/accordion_image_0001.surf";

    // SURF抽出用に画像をグレースケールで読み込む
    IplImage* grayImage = cvLoadImage(imageFile, CV_LOAD_IMAGE_GRAYSCALE);
    if (!grayImage) {
        cerr << "cannot find image file: " << imageFile << endl;
        return -1;
    }

    // キーポイント描画用にカラーでも読み込む
    IplImage* colorImage = cvLoadImage(imageFile, CV_LOAD_IMAGE_COLOR);
    if (!colorImage) {
        cerr << "cannot find image file: " << imageFile << endl;
        return -1;
    }

    CvMemStorage* storage = cvCreateMemStorage(0);
    CvSeq* imageKeypoints = 0;
    CvSeq* imageDescriptors = 0;
    CvSURFParams params = cvSURFParams(500, 1);

    // 画像からSURFを取得
    cvExtractSURF(grayImage, 0, &imageKeypoints, &imageDescriptors, storage, params);
    cout << "Image Descriptors: " << imageDescriptors->total << endl;

    // SURFをファイルに出力
    writeSURF(surfFile, imageKeypoints, imageDescriptors);

    // 画像にキーポイントを描画
    for (int i = 0; i < imageKeypoints->total; i++) {
        CvSURFPoint* point = (CvSURFPoint*)cvGetSeqElem(imageKeypoints, i);
        CvPoint center;  // キーポイントの中心座標
        int radius;      // キーポイントの半径
        center.x = cvRound(point->pt.x);
        center.y = cvRound(point->pt.y);
        radius = cvRound(point->size * 1.2 / 9.0 * 2.0);
        cvCircle(colorImage, center, radius, cvScalar(0,255,255), 1, 8, 0);
    }

    cvNamedWindow("SURF");
    cvShowImage("SURF", colorImage);
    cvWaitKey(0);

    // 後始末
    cvReleaseImage(&grayImage);
    cvReleaseImage(&colorImage);
    cvClearSeq(imageKeypoints);
    cvClearSeq(imageDescriptors);
    cvReleaseMemStorage(&storage);
    cvDestroyAllWindows();

    return 0;
}

たくさん後始末しなくちゃいけないので大変ですね・・・C++はこれだから嫌い。もし、何か解放し忘れてるのに気が付いたら教えてください。当初は、cvClearSeq()するのを知らず書いてなかったことがありました・・・

SURFを格納するファイルはSIFT特徴量の抽出（2009/10/24）で使ったSIFTファイルと同じ形式です。1行に1つのキーポイントの情報をずらずら出力してます。SURFの理論が詳しくわかってないのですが、サイズはSIFTのスケールのことかな？あとラプラシアンもあまりわかってないですが、1か-1の値を取ります。あとでキーポイントのマッチングを高速化するときに使うのでとりあえず保存しときます。

408 128   <-- キーポイントの数 SURFの次元数
138.821 11.208 14 1 0.00771829 0.00876939 ...     <-- 1つめのキーポイントのX座標, Y座標, サイズ, ラプラシアン, 以下128次元の特徴ベクトル
121.597 17.7526 17 -1 -0.00576097 0.00586741 ...  <-- 2つめのキーポイントの情報
...

SIFTと同様に抽出対象となる画像はグレースケール（モノクロ）で読み込みます。他の関数の使い方は上のサンプルで大体わかると思いますが、キーポイントの情報（imageKeypoints）とキーポイントの特徴ベクトルの情報（imageDescriptors）は別々のCvSeqに格納される仕様になってるようです。ただし、添え字は対応しています。たとえば、imageKeypointsのi番目のキーポイントの特徴ベクトルはimageDescriptorsのi番目に格納されてます。キーポイント情報はCvSURFPointという構造体にまとめられていて、特徴ベクトルはfloatの配列がそのまま入っているようです。SURFのパラメータは、cvSURFParamsで指定します。OpenCVのマニュアル（docフォルダにあるopencv.pdf）によると第1引数は特徴量の次元で0だと64次元、1だと128次元です。次元が大きいほうが正確なのかな？第2引数は閾値でキーポイントのhessianがこの値より大きい点だけ抽出されます。適切な閾値は300-500とのこと。大きくすると抽出される点が少なくなるのがわかります。

いくつかの画像で試してみます。SIFTと比べると背景と物体の間の線（エッジ）が抽出されやすいのかな？陰陽を見ると違いがよくわかります。エッジをキーポイントとするのはよくないそうなのでSIFTに比べて少し精度が落ちるってことでしょうか？もし詳しい人がいたらアドバイスください。

動画からSURFをリアルタイム抽出する

次は、上のプログラムを少し拡張して、Webカメラからキャプチャした動画からSURFをリアルタイムに抽出してみます。このプログラムは、Webカメラを接続してないと動きません。

#include <cv.h>
#include <highgui.h>
#include <iostream>

using namespace std;

int main(int argc, char** argv) {
    CvCapture* capture;

    // カメラを初期化
    if ((capture = cvCreateCameraCapture(0)) == NULL) {
        cerr << "cannot find camera" << endl;
        return -1;
    }

    // ウィンドウを生成
    cvNamedWindow("SURF");

    IplImage* captureImage = cvQueryFrame(capture);
    while (true) {
        CvMemStorage* storage = cvCreateMemStorage(0);
        CvSeq* imageKeypoints = 0;
        CvSeq* imageDescriptors = 0;
        CvSURFParams params = cvSURFParams(500, 1);

        captureImage = cvQueryFrame(capture);

        // グレースケールに変換
        IplImage* grayImage = cvCreateImage(cvGetSize(captureImage), 8, 1);
        cvCvtColor(captureImage, grayImage, CV_BGR2GRAY);

        // フレーム画像からSURFを取得
        cvExtractSURF(grayImage, 0, &imageKeypoints, &imageDescriptors, storage, params);

        // 画像にキーポイントを描画
        for (int i = 0; i < imageKeypoints->total; i++) {
            CvSURFPoint* point = (CvSURFPoint*)cvGetSeqElem(imageKeypoints, i);  // i番目のキーポイント
            CvPoint center;  // キーポイントの中心座標
            center.x = cvRound(point->pt.x);
            center.y = cvRound(point->pt.y);
            cvCircle(captureImage, center, 2, cvScalar(0,255,255), CV_FILLED);
        }
        cvShowImage("SURF", captureImage);

        // ループ内で作成したオブジェクトは始末
        cvReleaseImage(&grayImage);
        cvClearSeq(imageKeypoints);
        cvClearSeq(imageDescriptors);
        cvReleaseMemStorage(&storage);

        // ESCキーが押されたらループを抜ける
        int key = cvWaitKey(30);
        if (key == 27) {
            break;
        }
    }

    // 後始末
    cvReleaseCapture(&capture);
    cvDestroyAllWindows();

    return 0;
}

ループ内で生成したオブジェクトはループ内で後始末します。これ忘れるとものすごい勢いでメモリリークしてメモリを食いつぶしそのうちクラッシュします。これまたC++に慣れてないとときどきやっちゃいます。実行すると下の動画みたいにリアルタイムで抽出できます。静止画ではなく動画を対象にするとSIFTに比べてSURFが高速なのが実感できます。これからは特に注意書きしない限り、特徴点の抽出にはSURFを使っていきます。OpenCVに関数あるから楽だし（ボソッ）。

3日で作る高速特定物体認識システム (4) 特徴点のマッチング（2009/11/2）へ続きます。

参考文献

• SURF: Speeded Up Robust Features - SURFのホームページ、プログラムもあり
• Bay, H., Ess, A., Tuytelaars, T. and Gool, L. V.: SURF: Speeded Up Robust Features (PDF), Computer Vision and Image Understanding (CVIU), 110(3), pp.346-359 (2008) - SURFの原論文
• surfmex - SURFのMatlab実装
• SURF - Wikipediaの記事

このシリーズの記事一覧

• 3日で作る高速特定物体認識システム (1) 物体認識とは（2009/10/18）
• 3日で作る高速特定物体認識システム (2) SIFT特徴量の抽出（2009/10/24）
• 3日で作る高速特定物体認識システム (3) SURFの抽出（2009/10/30）
• 3日で作る高速特定物体認識システム (4) 特徴点のマッチング（2009/11/2）
• 3日で作る高速特定物体認識システム (5) 物体モデルデータベースの作成（2009/11/14）
• 3日で作る高速特定物体認識システム (6) 線形探索を用いた特定物体認識（2009/11/22）
• 3日で作る高速特定物体認識システム (7) 最近傍探索の高速化（2009/12/12）

■[OpenCV] 3日で作る高速特定物体認識システム (2) SIFT特徴量の抽出 20:39

3日で作る高速特定物体認識システム (1) 物体認識とは（2009/10/18）の続きです。

今回は、画像からSIFT (Scale-Invariant Feature Transform) という局所特徴量を抽出するところを作ってみようと思います。

（注）SIFTは特許がとられているので製品には使えません

SIFT特徴量の抽出

まずは、局所特徴量の代表ともいえるSIFTを試してみます。OpenCVにはSIFTを抽出する関数がなかったのでRob Hess氏がC言語で実装したライブラリを試してみます。内部でOpenCVを使っているので事前にOpenCVのインストールが必要です。実装はOpenCV 1.1でされているみたいですが、2.0でもちょっと手直しすると動きました。Rob Hess氏のホームページからSIFT Feature Detectorのzip版を落とします。解凍するといろいろファイルができますが、SIFTを抽出するのに最低限必要なのは、

  siftfeat.c - メイン関数
  sift.c sift.h - SIFT抽出用関数
  imgfeatures.c imgfeatures.h - SIFT特徴量を描画、保存する関数
  utils.c utils.h - ユーティリティ関数
  beaver.png - テスト画像

です。これらを自分のプロジェクトに追加します。OpenCV2.0でコンパイルするといくつかエラーが出た*1ので直します。imgfeatures.cの373行目を

cvEigenVV( &M, &V, &E, DBL_EPSILON, 0, 0 );

siftfeat.cの28, 29行目を

char* img_file_name = "beaver.png";
char* out_file_name  = "beaver.sift";

同じく72行目を

cvSaveImage( out_img_name, img, 0);

に修正します。まだいろいろ警告は出ますがとりあえずコンパイルできます。beaver.pngという画像で実行すると下のような画像が表示されます。

矢印の始点はキーポイントの位置、矢印の長さはスケール、矢印の向きは勾配強度がもっとも強い向きを表しているようです。SIFTの表示方法はキーポイントの点だけ表示、矢印で表示、円で表示などいろいろあるみたいですね。スケールや勾配強度が何を意味しているかはもう少し調査が必要です。当面は、SIFTのキーポイントの位置とそれぞれの点が128次元の特徴ベクトルで表せることだけおさえとけばいいかな？

それぞれのキーポイントの128次元ベクトルは、beaver.siftというファイルに保存されています。矢印1本につき128個の数値（128次元ベクトル）です。下のようなファイルです。

114 128                                                    <--- キーポイントの数と特徴ベクトルの次元数
101.350182 136.130755 40.169879 0.771085                   <--- 1つめキーポイントのX座標, Y座標, スケール, 向き
 0 0 0 0 3 1 0 0 2 23 46 15 18 3 0 0 6 20 13 1             <--- インデントされている数字が128次元ベクトル（20個で改行が入る）
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 88 36 0 0
 81 95 57 47 185 114 2 7 185 155 19 6 19 6 1 22 22 0 0 0
 0 0 0 1 0 0 0 0 37 8 0 0 91 12 0 1 185 144 11 35
 185 50 0 0 23 28 8 95 40 1 0 0 0 0 0 4 0 0 0 0
 0 0 0 0 11 5 0 0 4 2 0 0 49 20 0 0 1 0 0 1
 0 0 0 0 0 0 0 0
127.871536 71.100579 15.768590 -2.024589                   <--- 2つめのキーポイント
 1 2 2 72 63 12 1 1 133 93 1 4 2 7 4 44 133 115 0 0
 0 0 0 20 9 4 0 0 0 0 0 0 23 0 1 9 107 20 1 8
 133 5 0 0 0 1 5 133 132 14 0 0 0 0 8 133 14 1 0 0
 0 0 0 8 26 0 0 0 126 37 8 22 133 47 0 0 0 0 3 52
 131 41 0 0 0 0 2 36 1 0 0 0 0 0 0 2 2 0 0 0
 34 105 80 24 111 15 0 0 0 1 55 66 79 21 0 0 0 0 0 5
 0 0 0 0 0 0 0 0

上の画像の場合、128次元ベクトル114本によって表されるってことですね。この128次元ベクトル114本を物体モデルデータベースに登録しておくことになります。けっこう膨大な数値です。このファイルをよく観察するとどうやらキーポイントはスケールの大きな順（矢印が長い順）に並んでいることがわかります。これまたまだよくわかってないのですが、スケールが大きいキーポイントほど重要と解釈してよいのでしょうか？もし詳しい方がいたらアドバイスください。

キーポイントを点で描画

ここでちょっと改造してみます。キーポイントを矢印ではなく点表示にし、SIFTファイルが1行で1個のキーポイント情報を表すようにします。imgfeature.cの595行目あたりを

//	cvLine( img, start, end, color, 1, 8, 0 );
//	cvLine( img, end, h1, color, 1, 8, 0 );
//	cvLine( img, end, h2, color, 1, 8, 0 );

	cvCircle(img, cvPoint(start_x, start_y), 2, cvScalar(0,255,255,255), -1, 0, 0);

とし、同じファイルの523行目あたり（数字20個おきに改行を入れる部分）をコメントアウトします。

//			if( j % 20 == 0 )
//				fprintf( file, "\n" );

こうすると下のようにキーポイントが点になります。beaver.siftも1つのキーポイントの情報が1行にまとまります。スクリプトで処理するときはこっちの方が楽だと思います。

いくつかの画像でSIFTキーポイントを描画してみます。SIFTが画像のどういうところをキーポイントとしてるか大体想像つきますね。陰陽マークみたいにあまり模様がない画像ではキーポイントが少なくなりがちみたいです。

SIFTの詳細

SIFTのキーポイントは、画像の特徴をよく表している場所ってことはわかりました。では、どういう根拠に基づいてその点が選ばれたのか？128次元ベクトルがどういう風に求められたのか？はアルゴリズムの詳細を詳しく調べる必要があります。SIFTは結構複雑なことをやっていて、ガウシアン平滑化フィルタ、Difference of Gaussian (DoG)、局値検出、画像の微分、強度勾配など画像処理の理論がわかっていないと理解が難しい・・・私もまだおぼろげに理屈がわかってきた程度です。1つ思い違いをしてたのは、SIFTは色情報を使ってないってことかな。実際、グレースケール画像から抽出してます。いくつか参考にした文献を挙げておきます。これらを参考にそのうちC++版のSIFTを作ってみようと思います。

3日で作る高速特定物体認識システム (3) SURFの抽出（2009/10/30）へつづきます。

• SIFTチュートリアル - 藤吉先生によるSIFTの解説資料、図が豊富でわかりやすい。
• SIFTを1から組むブログ - 手順を追って少しずつ実装しているのでわかりやすい。途中でららら♪ってごまかしてるけど（笑）
• Rob Hess氏のSIFT特徴量抽出モジュールつかってみる - Rob Hess氏のプログラムの実行方法です。
• SIFT++ - 別の実装。こっちも試してみたい。
• SIFT for Matlab - Matlabによる実装。
• Scale-invariant feature transform - Wikipediaの記事。
• Keypoint detector - SIFTの提案者Lowe氏のデモ。
• Lowe, D. G.: Object recognition from local scale-invariant features (PDF), Proc. of IEEE Intl. Conf. on Computer Vision (ICCV), pp.1150-1157 (1999) - SIFTの原論文
• Lowe, D: Distinctive imagefeatures from scale-invariant keypoints (PDF), Proc. of Intl. Journal of Computer Vision (IJCV), 60(2), pp.91-110 (2004) - SIFTの原論文

このシリーズの記事一覧

• 3日で作る高速特定物体認識システム (1) 物体認識とは（2009/10/18）
• 3日で作る高速特定物体認識システム (2) SIFT特徴量の抽出（2009/10/24）
• 3日で作る高速特定物体認識システム (3) SURFの抽出（2009/10/30）
• 3日で作る高速特定物体認識システム (4) 特徴点のマッチング（2009/11/2）
• 3日で作る高速特定物体認識システム (5) 物体モデルデータベースの作成（2009/11/14）
• 3日で作る高速特定物体認識システム (6) 線形探索を用いた特定物体認識（2009/11/22）
• 3日で作る高速特定物体認識システム (7) 最近傍探索の高速化（2009/12/12）

■[OpenCV] 3日で作る高速特定物体認識システム (1) 物体認識とは 19:45

情報処理学会の学会誌『情報処理』の2008年9月号（Vol.49, No.9）に「3日で作る高速特定物体認識システム」という特集記事があります。OpenCVを用いた面白そうなプロジェクトなのでレポートにまとめてみようと思います。3日でできるかはわからないけど。

残念ながらこの記事はPDFを無料でダウンロードすることができません。なので会員以外で元記事が読みたい人は図書館でコピーする必要があるかも・・・また、2009年9月号の人工知能学会誌にも物体認識の解説「セマンティックギャップを超えて―画像・映像の内容理解に向けてー」があります。こちらも非常に参考になりますが同様にPDFが手に入りません・・・*1。他にもいくつかわかりやすい総説論文へのリンクを参考文献にあげておきます。

物体認識とは

物体認識（object recognition）は、画像に写っているものが何であるかを言い当てる処理で、特定物体認識と一般物体認識に分類されています。特定物体認識は、ある特定の物体と同一の物体が画像中に存在するかを言い当てる（identification）処理で、一般物体認識は、椅子、自動車、虎など一般的な物体のカテゴリを言い当てる（classification）処理です。特定物体認識はコンピュータの進歩で比較的容易になってきており商用利用もされているそうですが、本当に難しいのは一般物体認識で、近年の研究はこちらが中心でブームだそうです。

一般物体認識の難しさはセマンティックギャップと言われています。たとえば、同じ自動車というカテゴリでもいろんな形・色・種類の自動車がありますからね。人間はどんな色・形でもなぜか自動車とずばり言い当てられますが、コンピュータには自動車の本質的な意味がわからないため難しいわけです。というか人工知能の困難さは機械が意味（言葉でも画像でも音声でも）を理解できないことからきていると思います。近年はWeb上のタグ付き大量画像データと機械学習の手法を駆使して意味を自動学習させようとする研究が多い印象を受けています。

今回は表題の記事を参考に実装が比較的容易な特定物体認識を作ってみます。

物体認識システムの構成

各論に入る前にシステムの全体像を示しときます。

※http://www.m.cs.osakafu-u.ac.jp/IPSJ_3days/から引用

特定物体認識では、認識させたい物体をWebカメラなどで撮影して局所特徴量を抽出し、物体モデルデータベースにあらかじめ登録しておきます。検索は、クエリとして認識させたい画像を写して局所特徴量を抽出し、物体モデルデータベース中の各物体の局所特徴量との間で類似度を計算します。その中から類似度がもっとも高いものを認識結果として返すというのが全体の流れです。

表題の記事では、下の順でプロジェクトを進めようって感じです。

(1) 低速特定物体認識システムの構築（1日目）

局所特徴量の代表であるSIFT（Scale-invariant feature transform）特徴量を画像から抽出するモジュールを作成し、物体データベースモデルの検索にもっとも単純な線形探索で最近傍（類似度がもっとも近い点）を探す。線形探索は遅すぎて実用的ではないが、とりあえず動くシステムを作る。

(2) 近似最近傍探索を用いた高速化（2日目）

線形探索は遅すぎるので、物体モデルデータベースの局所特徴量をkd-treeやLocality Sensitive Hashing（LSH）でインデキシングして最近傍探索を高速化する。

(3) ユーザインタフェースの構築（3日目）

物体モデルデータベースの作成や認識結果の表示を行うインタフェースを作成する。

3日で進めるのはけっこうハードです・・・

局所特徴量

特定物体認識システムでは、データベースに登録されたのと同一の物体のみ認識させることができます。面白いのはまったく同じ画像じゃなくても認識できることですかね。物体を撮影する向き、サイズが多少変わっていたり、一部が隠れていてもきちんと認識できます。これは、画像から抽出した多数の局所特徴量を使うことで可能となります。画像の局所特徴量によるモデル化は、この分野を大きく発展させたブレークスルーだそうです。下は局所特徴量の代表であるSIFTの例です。

左は元のレナさん画像です。右はサイズを80%に縮小し（スケール）、90度回転させ、明度を落とし（照明）、周囲を隠蔽（オクルージョン）した画像です。個々の黄色い点がキーポイント（Interest point、特徴点など呼び方はいろいろ）で、各キーポイントがそれぞれ128次元のSIFT特徴量を持ちます。この画像だと大体1000個くらいのキーポイントが得られます。比べてみるとわかるように画像が多少変化していてもキーポイントの位置は大体同じでSIFT特徴量の値も近い（後で出力します）ことがわかります。SIFTは画像のスケール変化、照明変化、平行移動、回転、隠蔽（オクルージョン）に対して頑健といわれる理由です。つまり、左のレナさんの局所特徴量を物体モデルデータベースに登録しておけば、多少変化してしまった右のレナさんも似た局所特徴量を持つためちゃんと認識できるってわけですね。

類似画像検索システムを作ろう（2009/10/3）では、カラーヒストグラムという特徴量を使いましたが、これは大域（global）特徴量といって画像全体から1つの特徴量を得ます。一方、局所（local）特徴量は画像全体からではなく、画像の局所から多数（上の例では約1000個）の特徴量を得ています。これはけっこうすごい考え方の転換だと思います。

今日は疲れたのでここまで。3日で作る高速特定物体認識システム (2) SIFT特徴量の抽出（2009/10/24）につづきます。

参考文献

• 黄瀬浩一, 岩村雅一: 3日で作る高速特定物体認識システム, 情報処理, Vol.49, No.9, pp.1082-1089 (2008)
• 井出一郎, 柳井啓司: セマンティックギャップを超えて―画像・映像の内容理解に向けて―, 人工知能学会誌, Vol.24, No.5, pp.691-699 (2009)
• 柳井啓司: 一般物体認識の現状と今後（PDF）, 情報処理学会研究報告CVIM, Vol.93, pp.121-134 (2006)
• 3日で作る高速特定物体認識システム - プロジェクトのサイト
• 意味的マルチメディア処理特論 - 講義資料など多数
• コンピュータビジョン特論 - 講義資料

このシリーズの記事一覧

• 3日で作る高速特定物体認識システム (1) 物体認識とは（2009/10/18）
• 3日で作る高速特定物体認識システム (2) SIFT特徴量の抽出（2009/10/24）
• 3日で作る高速特定物体認識システム (3) SURFの抽出（2009/10/30）
• 3日で作る高速特定物体認識システム (4) 特徴点のマッチング（2009/11/2）
• 3日で作る高速特定物体認識システム (5) 物体モデルデータベースの作成（2009/11/14）
• 3日で作る高速特定物体認識システム (6) 線形探索を用いた特定物体認識（2009/11/22）
• 3日で作る高速特定物体認識システム (7) 最近傍探索の高速化（2009/12/12）

■[OpenCV] EclipseでOpenCV 21:35

OpenCVのプログラムを統合開発環境EclipseとC/C++コンパイラMinGWで開発するためのセットアップ方法をまとめておきます。Windows環境でVisual StudioをインストールしなくてもOpenCVプログラムを開発できます。

1. MinGWのインストール

ここからMinGW-5.1.6.exe（2009/10/15現在）をダウンロード。インストール途中でg++ compilerとMinGW Makeオプションをチェック。

2. Eclipseのインストール

ここからEclipse IDE for C/C++ Developers（Eclipse CDT）をダウンロード。解凍するとeclipseフォルダができるのでC:\eclipseなどに移動。

3. OpenCVのインストール

ここからOpenCV-2.0.0a-win32.exe（2009/10/15現在）をダウンロード。デフォルトでは、C:\OpenCV2.0にインストールされます。

4. 環境変数の設定

環境変数のPATHにC:\MinGW\bin;C:\OpenCV2.0\bin;を追加。OpenCVのbinを追加しないとdllがリンクできずOpenCVプログラムがクラッシュするので注意。

5. 新規プロジェクトを作成

Eclipseで新規C++プロジェクトを作成。プロジェクト名（ここではtest）を入力。Project typeでEmpty Project、ToolchainsでMinGW GCCを選択。

6. includeとlibフォルダの設定

コンパイラにOpenCVのヘッダファイル（include）とライブラリファイル（lib）の場所を教えます。Project Explorerのtestプロジェクトを右クリックしてPropertiesにてC/C++ General > Paths and Symbolsを選択。Configurationで[All configurations]を選択。Includes > C++ Source Fileを選択してC:\OpenCV2.0\include\opencvをAdd。Library Pathsを選択してC:\OpenCV2.0\libをAdd。

（注）CとC++の設定は別々なので注意。もし、拡張子が.cのCプログラムを追加するならIncludes > C Source Fileにも同様にC:\OpenCV2.0\include\opencvを追加する必要があります。

7. リンクするライブラリの設定

C/C++ Build > Settingsを選択。Configurationで[All configurations]を選択。MinGW C++ Linker > LibrariesのLibrariesにcv200、highgui200、cxcore200、ml200をそれぞれ追加。

（注）CとC++の設定は別々なので注意。もし、拡張子が.cのCプログラムを追加するならC Source Fileにも追加する必要がある。

（注）5から7の設定はプロジェクトを作成するたびに行います。面倒ですが・・・

8. ソースファイルを作成

プロジェクトにソースファイルを追加してコードを書きます。たとえば、main.cpp（C++プログラムです）を追加して以下のコードを入力します。Eclipseは入力補完がすごく充実してます。たとえば、cvと打ち込んだ状態でCtrl+Spaceを押すと関数候補が引数の説明付きでずらっと並ぶので↑↓矢印で選択してEnterを押すとかってに入力してくれます。

#include <cv.h>
#include <highgui.h>

int main(int argc, char **argv) {
    // 画像をロード
    IplImage *img = cvLoadImage("C:/OpenCV2.0/samples/c/lena.jpg", CV_LOAD_IMAGE_COLOR);

    // ウィンドウを作成
    cvNamedWindow("image", CV_WINDOW_AUTOSIZE);
    // ウィンドウに画像を表示
    cvShowImage("image", img);

    // キー入力待機
    cvWaitKey(0);

    // 後処理
    cvReleaseImage(&img);
    cvDestroyWindow("image");

    return 0;
}

9. ビルド

メニューのProject > Build AllかCtrl+Bでビルド。成功すればProject ExplorerにBinaries/test.exeができます。もし失敗したらソースコード上に警告やエラー位置と内容が表示されるので修正すべき箇所が一目でわかります。

（注）__exchange_and_addに関連したコンパイルエラーが出る場合は以下の修正を加えます（OpenCV InstallGuideより）。どうやらVisual Studioが入っていない環境で出るエラーみたいです。

(1) C:\OpenCV2.0\include\opencv\cxoperations.hppをエディタで開く。改行コードの関係でメモ帳では崩れます。Eclipseで開いてください。

(2) 66から68行目（OpenCV2.0.0の場合）あたりの

  #else
    #include <bits/atomicity.h>
    #if __GNUC__ >= 4

の部分を

  #else
    #include <bits/atomicity.h>
    #if __GNUC__ >= 4 || __MINGW__

に修正する。__MINGW__はMINGWの前と後にアンダースコアが2つずつあります。

(3) Eclipseでtestプロジェクトを右クリックしてPropertiesを開きます。C/C++ General > Paths and Symbols > Symbolsを開いてConfigurationで[All Configuration]、LanguagesでC++ Source Fileを選択。Nameに__MINGW__、Valueは空のエントリを追加。

10. 実行

Ctrl+F11でプログラムを実行です。レナさんが表示されれば大成功！