読者です 読者をやめる 読者になる 読者になる

MyEnigma

とあるエンジニアのブログです。#Robotics #Programing #C++ #Python #MATLAB #Vim #Mathematics #Book #Movie #Traveling #Mac #iPhone

高速最近傍点探索ライブラリFLANNを使ってみる

はじめてのパターン認識

はじめてのパターン認識

目次

はじめに

あるデータ群の中で、

最も近いデータを選ぶ

最近傍探索は、

様々なアルゴリズムで使用されますが、

それを高速に実現するには、

データ構造から工夫する必要があります。

 

今回は、上記の目的を達成するための

C++ライブラリFLANNを紹介したいと思います。

 

最近傍点探索とは

最近傍点探索とは、

N次元の空間内の、M個のデータの中から、

ある一つのデータ(クエリデータ)に

最も近いデータを計算することです。

 

この最近傍点探索は、

空間の次元Nが増えたり、

データ数Mが増えると

非常に計算時間がかかってしまうので、

この探索を高速化する様々な方法が研究されてきました。

 

FLANN (Fast Library for Approximate Nearest Neighbors)

FLANN(Fast Library for Approximate Nearest Neighbors)は、

C++で書かれた高速最近傍点探索ライブラリです。

下記のGithubリポジトリでソースコードが公開されています。

github.com

 

作者の研究によると、

これまで開発された様々な最近傍点探索のアルゴリズムと比べて、

約10倍高速とのことです。

また、このFLANNはC++で書かれていますが、

簡単にCやMATLAB, Pythonで利用することができます。

 

加えて、FLANNはPoint Cloud Library(PCL)や、

Open CVなどでも使用されているようです。

 

また、FLANNはLinux上で動くことはテストされているようですが、

ユーザマニュアルによると、WindowsやMacでも利用できるようです。

ライセンスはBSDライセンスなので、

商用利用もしやすいと思います。

 

今回の記事では、C++とPythonで

このFLANNを利用する方法をまとめたいと思います。

 

インストール方法

下記のリンクのGetting FLANNから

ソースコードのZipをDLしましょう。

 

続いて、下記のようにbuildディレクトリを作って、

cmakeでビルドします。

$ cd flann-x.y.z-src

$ mkdir build

$ cd build

$ cmake ..

$ make

 

以上でビルドは完了です。

  

C++でFLANNを使う

自分の環境では(Ubuntu14.04)、

ユーザマニュアルにある通りに、

コンパイルしようとした所、

エラーでコンパイルすることができませんでした。

flann_manual-1.8.4.pdf

  

下記のように、

flannのヘッダが/usr/includeにダウンロードされていたことと、

hdf5というライブラリが必要だったようなので、

下記のようなコマンドでコンパイルすることができました。

$ g++ ${cwd}/sample1.cpp -I /usr/include/ -lhdf5 -o test

 

下記のGithubリポジトリに、

コンパイル用シェルスクリプトと、

サンプルコードを公開しました。

github.com

こちらのコードを元に、

FLANNの基本的な使い方について説明しようと思います。

 

Radius Searchによる最近傍点探索

下記のコードは、

Radius Searchによる最近傍点探索のサンプルコードです。

総当り法と計算時間の比較をしています。

 

このシミュレーションは、

3次元のデータを10000個作成し、

その中からある指定したデータqueryに

最も距離が近いデータを選びます。

 

自分の環境の場合、

総当りは0.8秒

FLAANのkd-treeのインデックス探索では0.2秒と、

約計算時間が1/4になりました。

非常に高速ですね。

 

/**
 *  radius search sample
 * */
#include <flann/flann.hpp>
#include <flann/io/hdf5.h>

#include <stdio.h>
#include <random>
#include <chrono>
#include <math.h>

using namespace flann;
using namespace std;

void ShowMatrix(const Matrix<int> &mat){
  for(int i=0;i<mat.rows;i++){
    for(int j=0;j<mat.cols;j++){
      cout<<mat[i][j]<<",";
    }
    cout<<endl;
  }
}


void ShowMatrix(const Matrix<float> &mat){
  for(int i=0;i<mat.rows;i++){
    for(int j=0;j<mat.cols;j++){
      cout<<mat[i][j]<<",";
    }
    cout<<endl;
  }
}

int BruteForceSearch(
    const Matrix<float> &dataset,
    const Matrix<float> &query
    ){
  int ind=-1;
  double dist=10000.0;

  for(int i=0;i<dataset.rows;i++){
      // cout<<mat[i][j]<<",";
      double dx=query[0][0]-dataset[i][0];
      double dy=query[0][1]-dataset[i][1];
      double dz=query[0][2]-dataset[i][2];

      double tdist=sqrt(dx*dx+dy*dy+dz*dz);

      if(tdist<dist){
        dist=tdist;
        ind=i;
      }
  }

  cout<<"Brute dist:"<<dist<<endl;
  return ind;
}

int main(int argc, char** argv)
{
    cout<<"start"<<endl;
    int nDim = 3;
    int nData = 10000;

    Matrix<float> dataset(new float[nDim*nData], nData, nDim);

    random_device rd;
    mt19937 mt(rd());
    
    //datasetを作成
    uniform_real_distribution<float> score(-10.0,10.0);
    for(int i=0;i<nData;i++){
      for(int j=0;j<nDim;j++){
        dataset[i][j]=score(mt);
        // cout<<dataset[i][j]<<endl;
      }
    }

    Matrix<float> query(new float[nDim*1],1, nDim);
    query[0][0]=10.0;
    query[0][1]=2.0;
    query[0][2]=10.0;
    // cout<<"query:";
    // ShowMatrix(query);

    Matrix<int> indices(new int[query.rows],query.rows,1);
    Matrix<float> dists(new float[query.rows],query.rows, 1);

    // cout<<query.rows<<endl;
    Index<L2<float> > index(dataset, flann::KDTreeIndexParams(4));
    index.buildIndex();                                                                                               
    chrono::system_clock::time_point  start, end; 

    //Bruteforce
    start = chrono::system_clock::now();
    int ind=BruteForceSearch(dataset,query);
    cout<<"brute:ind:"<<ind<<endl;
    end = chrono::system_clock::now();
    double msec = (double)chrono::duration_cast<chrono::nanoseconds>(end-start).count()/1000.0;
    cout<<"time:"<<msec<<endl;

    start = chrono::system_clock::now();
    cout<<index.radiusSearch(query, indices, dists, 100.0,flann::SearchParams(128))<<endl;
    end = chrono::system_clock::now();
    msec = (double)chrono::duration_cast<chrono::nanoseconds>(end-start).count()/1000.0;
    cout<<"time:"<<msec<<endl;

    cout<<"indices:";
    ShowMatrix(indices);
    cout<<"dists:";
    ShowMatrix(dists);
    for(int i=0;i<indices.cols;i++){
      // cout<<indices[0][i]<<endl;
      cout<<dataset[indices[0][i]][0]<<",";
      cout<<dataset[indices[0][i]][1]<<",";
      cout<<dataset[indices[0][i]][2]<<endl;
    }

    delete[] dataset.ptr();
    delete[] query.ptr();
    delete[] indices.ptr();
    delete[] dists.ptr();
    
    return 0;
}

参考資料

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

はじめてのパターン認識

はじめてのパターン認識