DEIM Forum 2017 E7-3

SIFT 特徴量を用いた画像検索の GPU による高速化
草村

優太†

天笠

俊之††

北川 博之††† 筑波大学大学院システム情報工学研究科 〒 305–8573 茨城県つくば市天王台 1–1–1†† 筑波大学計算科学研究センター 〒 305–8573 茨城県つくば市天王台 1–1–1
E-mail: †kusamura@kde.cs.tsukuba.ac.jp, ††{amagasa,kitagawa}@cs.tsukuba.ac.jp
あらまし

画像検索とは，画像集合中からクエリ画像と視覚的に類似した画像を見つけ出す処理であり，様々なアプ

リケーションで利用されている．画像間の類似度を計算するため，画像検索では画像から抽出される特徴ベクトルを用いる手法が主流である．しかし，画像検索では扱う特徴ベクトル数が膨大であり，特徴ベクトルが高次元であるため，多大な処理コストを要する．そのため，高速な画像検索の実現は非常に重要である．本研究では，GPU (graphics

processing unit) を用いた画像検索の高速手法を提案する．その際，データ量の削減と処理の高速化のために，提案手法では特徴ベクトルの圧縮を行う．その上で GPU 上での計算に適したデータ構造及びアルゴリズムを用いて画像検索を行う．さらに，画像データ及び動画データを用いた実験により，提案手法の性能を評価する．
キーワード

GPU，画像検索，SIFT 特徴量

1. は じ め に
画像データベース中からクエリ画像と視覚的に類似した画像を検索する処理のことを画像検索と呼ぶ [4]．画像データは普遍的に利用されるため，画像検索は様々なアプリケーションで利用されている．

可能である．しかし，画像データベースでは特徴ベクトル数が莫大となり，さらに，特徴ベクトルが高次元であるため，画像検索は多大な処理コストを要する．そのため，画像検索の高速化は非常に重要である．
一方，近年では莫大なデータを扱う処理の高速化のためにGPU (graphics processing unit) を用いた並列コンピューティ

画像検索を適用可能なアプリケーションとして，スマートングが注目されている．GPU は従来グラフィック処理に用いらフォン等を利用した歩行者の位置推定 [17] があげられる．このれるプロセッサであるが，多数のコアを搭載し高い並列処理性

アプリケーションは利用者として視覚障害者を想定し，利用者を持つという特徴がある．この特徴を活かして GPU を汎用計算の日常生活における移動を支援する．事前に普段よく移動するに用いる技術のことを GPGPU (general-purpose computing

経路を位置情報付きのビデオ映像で撮影しておき，各フレーム

on graphics processing units) と呼び，様々なアプリケーショとその位置情報をデータベースに格納しておく．位置推定は，ンにおいて処理の高速化に貢献している [14]．

歩行者が撮影した現在のカメラ映像に最も類似した画像をデー

そこで本稿では，SIFT 特徴量を用いた画像検索の GPU に

タベース上から検索し，その位置情報を取得することで行う．

よる高速化手法を提案する．提案手法では，特徴ベクトルの検索精度を向上させるためには，同じ経路であっても異なる気圧縮のために LSH (locality sensitive hashing) [11] を用いる．

象条件やライティングの条件が異なる複数の映像を格納しておLSH は類似したベクトルを高確率で同一のハッシュ値に変換すくことが望ましい反面，検索時には各映像フレーム毎に位置をるハッシュ法であり，類似検索の高速化等に利用されている．

特定する必要があるため，大量の画像に対して高速な検索を行提案手法は，データベース構築とクエリ処理の二つの処理にう必要がある．また，検索結果は最もマッチする上位 1 件のみ大別される．まず，データベース構築では，画像データベースが必要であり，Top-k 検索のような複雑な処理は不要である．

から抽出された特徴ベクトルを GPU での処理に適したデータ

このようなアプリケーションにおいて，膨大なビデオ映像と構造に変換し，GPU への転送を行う．また，クエリ処理は与位置情報が必要となるため，データベースはサーバー上で保持えられたクエリ画像を GPU 上で検索する処理である．

されることが想定される．位置推定の際には，まず歩行者が撮

また，本研究では提案手法の性能を評価するために実験を行っ影したカメラ映像をクエリ画像としてサーバーに送信する．そた．画像および動画データセットを用いた実験により，CPU をの後，サーバー上でデータベースとの照合を行い，得られた位用いた検索手法と比較して約 3.1 倍の高速化が可能であること

置情報をユーザーに送信する．このように，サーバー上でクエを示した．また，歩行者の位置推定アプリケーションへの適用リを処理することにより，高速な処理が可能となりリアルタイを想定した精度評価を行い，62 % の正解率を実現できること

ムに位置情報を取得可能となる．

を示した．

画像の特徴を表すベクトルのことを特徴ベクトルといい，画

本稿の構成を以下に示す．まず，2. 節にて，本研究に必要と像検索では画像間の類似度計算に特徴ベクトルを用いる手法がなる前提知識について説明する．次に，3. 節にて，本研究の関主流である．特徴ベクトルは数千程度のベクトルを用いて画像連研究を紹介する．その後，4. 節にて，本研究の提案手法につの特徴を表現し，画像検索はこれらを比較することにより実現いて説明し，その性能評価のための実験について 5. 節にて説明する．最後に，6. 節にて，まとめと今後の展望について説明

確率で同一のハッシュ値に変換するハッシュ法の一つである．

する．

画像検索において莫大な画像を扱う場合，特徴ベクトル数も莫2. 前 提 知 識

大となる．そのため，LSH は特徴ベクトル圧縮は近傍検索高速化に利用される [3, 8]．

本節では，本研究に必要となる前提知識について説明する．

LSH では，複数のハッシュ関数 hi () (0 <= i < k) からなる

2. 1 節にて，画像検索について説明し，2. 2 節にて GPU およハッシュ関数 H () = (h0 (), h1 (), . . . , hk −1 ()) を用いてベクトル

び GPU コンピューティングについて説明する．

をハッシュ値に変換する．つまり，ベクトル v は以下の式によ2. 1 画 像 検 索
画像データベース中からクエリ画像と視覚的に類似した画像を検索する処理のことを画像検索と呼ぶ [4]．画像検索は様々

り，サイズ k のハッシュ値 v ′ に変換される．

v ′ = H (v) = (h0 (v), h1 (v), . . . , hk−1 (v))

なアプリケーションで利用されているが，本研究では画像検索

ハッシュ関数 hi () (0 <= i < k) は，確率的に生成され，その生成を利用した歩行者の推定のアプリケーションを想定し，本研究

方法は圧縮するベクトル空間の距離指標によって異なる．現在

で扱う画像検索を以下のように定義する．まず，入力として画

では，lp ノルム，コサイン類似度，ハミング距離など，様々な像データベースと一つのクエリ画像が与えられる．その際，ク

距離指標に対するハッシュ関数が提案され利用されている [2]．

エリ画像と最も視覚的に類似した画像を画像データベース中から見つけ出し，その画像を結果として出力する．

SIFT 特徴ベクトル間の距離は l2 ノルムが用いられることが標準的であるため，本研究では l2 ノルムに対する LSH [9] を用

画像間の類似度を計算するために，近年の画像検索では局所いる．l2 ノルムに対する LSH では，三つの変数 a，b，W に基特徴量を用いる手法が主流である．特に，SIFT (scale-invariant

づいてハッシュ値を計算する．a は d (d = |v|) 次元のベクトル

feature transform) 特徴量 [12] は，最も基本的な特徴量としてであり，ガウス分布から選択された d 個の確率変数を要素とす

画像認識の分野で広く用いられている．そこで，本研究ではる．W は W > 0 となる任意の実数であり，b は閉区間 [0, W ]SIFT 特徴量を用いた画像検索の高速化手法を提案する．

の一様分布から選択される実数である．ハッシュ関数 hi () はこ

2. 1. 1 SIFT 特徴量

れらを用いて，以下の式で表される．

SIFT 特徴量 [12] は局所特徴量の一つであり，数千程度の128 次元ベクトルで画像の特徴を表現する．SIFT 特徴量は照明変化，拡大縮小，回転に対して強固であるという特徴を持つ．
そのため，視覚的に類似した画像を見つける処理である画像検索に適している．

2. 1. 2 SIFT 特徴量を用いた画像検索
SIFT 特徴量を用いた単純な画像検索 [1] は以下のように実現することができ，その処理はデータベース構築とクエリ処理の二つの処理に大別される．
まず，データベース構築について述べる．画像データベースが与えられると，画像データベース中の各画像から SIFT 特徴
量を抽出する．次に，抽出されたそれぞれの特徴ベクトル f と，
そのベクトルが抽出された画像 ID を示す値 id のタプル (f , id )を，データベース D に格納する．このとき，id (f ) を特徴ベク
トル f に対応する画像 ID として定義する．
また，クエリ処理では D に基づいて，与えられたクエリ画
像と最も類似した画像を返す．クエリ画像が与えられると，はじめにクエリ画像から特徴量 Q を抽出する．次に，クエリのそれぞれの特徴ベクトル q ∈ Q に対して，最も距離が近い D 中のベクトルに対応する画像 ID を得る．つまり，以下の式で表される画像 ID を得る．

id (fi ), s.t.∀fi , fj ∈ D ∧ dist(q, fi ) <= dist(q, fj )ただし，dist() は二つのベクトル間の距離である．この処理をクエリのそれぞれの特徴ベクトル q に対して行い，最も多く得られた画像 ID に対応する画像が検索結果となる．

2. 1. 3 LSH 及び SIFT 特徴量を用いた画像検索
LSH (locality sensitive hashing) [11] は，類似ベクトルを高hi (v) = ⌊

a·v+b
⌋
W

これは，元空間が a に直行する複数の超平面により等幅の部
分空間に分割され，同一部分空間に属しているベクトルが同一
ハッシュ値となることを意味する．この時，W の値はこのハッシュ関数のパラメータとなり，同一ハッシュ値に変換される部
分空間の領域サイズを決定する．具体的には，W が大きいほど領域サイズが大きくなり，同一ハッシュ値に変換されるベク
トル数が多くなる．そのため，LSH を近傍検索に用いた際のあいまい性が高くなる．

LSH 及び SIFT 特徴量を用いた画像検索 [16] は，2. 1. 2 節の手法を拡張することで実現できる．まず，データベース構築
では，特徴ベクトルと画像 ID のタプル (fi , id ) を得た後に各特徴ベクトルに LSH を適用する．すなわち，特徴ベクトル fiと LSH のハッシュ関数 H () に対して，ハッシュ値 fi′ = H (fi )を得る．その後，データベース D ′ にハッシュ値をキーとしそれに対応する画像 ID の多重集合を値とする連想配列を格納す
る．このとき，異なる特徴ベクトルから同一のハッシュ値が得られた場合は，それらに対応するすべての画像 ID を保持する
多重集合が連想配列の値となる．
クエリ処理では，クエリ画像から特徴量 Q が抽出されるとデータベース構築と同様に各特徴ベクトル q ∈ Q に LSH を適用して q ′ ∈ Q′ とする．ここで，得られた各ハッシュ値 q ′ に対して，D′ のキーから q ′ = f ′ となるような f ′ を検索し，対応する画像 ID を取得する．この処理で最頻出の画像 ID が検索結果となる．なお，LSH を用いているため，完全一致するハッ
シュ値の検索が特徴ベクトルの近傍検索の代用となる．

2. 2 GPU コンピューティング

を提案した．このフレームワークでは，特徴ベクトルを量子化

GPGPU (general-purpose computing on graphics process-

し複数の特徴ベクトルを同一の特徴ベクトルとみなす．この量

ing units) は，従来グラフィック処理向けに利用される GPU子化されたベクトルは visual words と呼ばれ，これらを用いて(graphics processing unit) を汎用計算に用いる技術のことで文書検索に基づいた方法でクエリ画像の検索を行う．

ある．GPU には，CPU と比較して多くのコアが搭載されていGPU を用いた画像検索の高速化に関する研究も行われているため，高い並列処理性を持つという特徴がある．そのため，る．Cevahir らの提案した手法 [6] は，2-ステップマッチング

GPGPU は科学技術計算やデータベースクエリ処理など様々なにより特徴ベクトルの検索を行う．2-ステップマッチングでは，

アプリケーションにおいて処理の高速化に貢献している [14]．

まず階層的 k-means を用いて特徴ベクトルをクラスタリング

GPU のアーキテクチャは独特であるため，不適切な実装やしておく．そしてクエリ画像が与えられると，クエリ画像の特アルゴリズムは CPU よりも悪いパフォーマンスを示す．その徴ベクトルに対して最も近いクラスタを検索し，さらにクラス

ため，GPU の特性を理解した上で適切なプログラムを組む必タ内で最近傍のベクトルを検索する．これらの処理を GPU 上要がある．本研究では，NVIDIA 社の GPU とその統合開発環

で行うことにより，画像検索の高速化を実現した．

境である CUDA

（注 1）

を用いる．

特徴ベクトルを扱うアプリケーションにおいては，その空間

GPU は一般的に，プロセッサとメモリにおいて階層的な構

コストや転送コストの削減が重要であり，SIFT 特徴ベクトル造を持つ．まず，GPU にはグローバルメモリとストリーミンの圧縮方法に関する研究も多く行われている．Chandrasekhar

グマルチプロセッサが搭載されている．グローバルメモリはらはハッシュ法，変換符号化，ベクトル量子化の三つに大別さGPU 上で最も大きなメモリであり，ストリーミングマルチプ

れる圧縮手法を比較し，その性能を評価した [7]．本研究で用ロセッサにはシェアードメモリとレジスタが搭載されている．

いる LSH はハッシュ法の一種であり，特徴ベクトルの圧縮にシェアードメモリは最大で 64 KB と比較的小さいメモリだが，有効である．

同一のストリーミングマルチプロセッサ内からの高速なアクセ

Carraher らは，SIFT 特徴ベクトルの高速な近傍検索手法でスが可能である．レジスタはシェアードメモリよりもさらに小ある LSH-NN を提案した [5]．LSH-NN では，LSH と GPU をさく高速なメモリであり，それぞれのスカラプロセッサからの用いることにより特徴ベクトルの近傍検索を高速に行う．しかみのアクセスが可能である．

し，LSH-NN で実現されたのは特徴ベクトルの近傍検索の高速

また，GPU はその処理モデルも階層的である．GPU には，
スレッド，ブロック，グリッドと呼ばれる三つの処理モデルが存在する．ブロックは数百程度のスレッド集合であり，グリッドは最大で数万個のブロックからなる．スレッド，ブロック，グリッドはそれぞれスカラプロセッサ，ストリーミングマルチプ
ロセッサ，デバイスで処理される．
本研究ではいくつかのデータ並列アルゴリズムプリミティブを用いる．データ並列アルゴリズムプリミティブは並列計算で化であり，特徴ベクトルを用いた画像検索の高速化は実現されていない．

4. 提 案 手 法
提案手法では，2. 1. 3 節で述べた LSH 及び SIFT 特徴量を用いた画像検索の高速化を実現する．その処理は，2. 1. 2 節及び 2. 1. 3 節と同様に，データベース構築とクエリ処理の二つの処理に分けられる．

頻出し，並列処理により高速な計算が可能である．また，GPU

本研究で高速化する処理は，クエリ処理の内，特徴ベクトル

向けのライブラリ（注 2） も存在し，利用することで高速な処理がを抽出してから画像 ID を返すまでの処理を対象とする．特に，可能である．本研究では，以下のデータ並列アルゴリズムプリ

LSH による特徴ベクトルの圧縮とデータベースとの照合処理のミティブを利用する．

高速化を図る．また，データベース構築では，データベースとreduce: 配列 [a0 , a1 , . . . , an−1 ] と二項演算子 ⊕ に対して，スの照合処理の高速化を実現するため，GPU 上での処理に適し

カラ値 a0 ⊕ a1 ⊕ . . . ⊕ an−1 を返す．

たデータベースに変換する．

scan: 配列 [a0 , a1 , . . . , an−1 ] と二項演算子 ⊕，単位元 I に対

4. 1 データベース構築

して，配列 [I , a0 , a0 ⊕ a1 , . . . , a0 ⊕ a1 ⊕ . . . ⊕ an−1 ] を返す．

データベース構築では，画像データベースを GPU に適した

sort: 配列 [a0 , a1 , . . . , an−1 ] を規則に基づいてソートする．

データ構造に変換した上で GPU に転送する．その処理では，これらのデータ並列アルゴリズムプリミティブは，全てライブ

まず 2. 1. 3 節と同様にデータベース D ′ を作成する．その後，ラリが利用可能である（注 3）．

D′ を元に以下の配列を作成し GPU に転送する．これは，配列

3. 関 連 研 究
Sivic らは BoVW (bag-of-visual-words) フレームワーク [15]が GPU 上での処理に適したデータ構造であるためである．

f val : データベース中の全ての fi′ を要素とする配列．要素は昇順でソートして格納される．

id val : fi′ 中の全ての画像 ID を f val に対応する順番で格納（注 1）：Parallel Programming and Computing Platform | CUDA |した配列．

NVIDIA www.nvidia.com/cuda

id ptr : id ptr [i] が id val の先頭レコードからのオフセット（注 2）：CUDA Toolkit Documentation http://docs.nvidia.com/cuda/（注 3）：Thrust :: CUDA Toolkit Documentation http://docs.nvidia.

を示す配列．

com/cuda/thrust/

クエリ処理で f val に対する二分探索を行うため，f val はデーAlgorithm 1

Algorithm 3

Input: D ′ ただし，(fi′ , idi ) ∈ D ′ ，idij ∈ idi

Input: Q ′ , f val, id ptr , id val

Output: f val, id val ，id ptr

Output: M

1: D′ を fi′ の値に基づき昇順ソート

1: for i = 0 : Q ′ .size() − 1 do in thread parallel

2: oﬀset = 02:

index = binarySearch(qi′ , f val)

3: id ptr [0] = 0

3:

if index == −1 then

4: for i = 0 : rowsof (D ′ ) − 1 do

4:

5:

f val[i] = fi′

6:

for j = 0 : sizeof (idi ) − 1 do

7:

id val[oﬀset + +] = idij
id ptr [i + 1 ] = oﬀset

8:

/

2

P[i] = −1

5:

else

6:

01

7:
8:
9:

N [i] = 0else
N [i] = id ptr [index + 1] − id ptr [index ]10:

Algorithm 2

11: F = scan(N )Input: Q ただし，qij ∈ Q

12: for i = 0 : Q ′ .size() − 1 do in block parallel′
Output: Q′ ただし，qij∈Q

13:

1: for i = 0 : Q.size() − 1 do in block parallelfor j = 0 : N [i] − 1 do in thread parallel

14:

for j = 0 : k − 1 do in thread parallel

2:

P[i] = id ptr [index ]if index == −1 thenM [F [i] + j] = id val[P[i] + j]′qij
= hk (qij )

3:

1

1

41 26
5317

0 ,

531 260 ,,,

図 1 データベース D ′ に基づく配列の作成．

タを昇順にソートする．また，ハッシュ値に対応する画像 ID

図2

配列 P，N，F の作成．

4. 2. 2 データベースとの照合の集合を二つの配列 id val および id ptr で表現する．これら

データベースとの照合では，クエリから得られたハッシュ値の配列の作成処理は，Algorithm 1 で表される．また，図 1 はと同一のハッシュ値を f val から検索し，一致したハッシュ値

配列作成処理の例である．

に対応する画像 ID を id val と id ptr を元に取得する．この

4. 2 クエリ処理

処理で，最も多く得られた画像 ID を検索結果とする．

クエリ処理では，データベース中からクエリ画像に最も類似

この処理は二分探索と atomicAdd を用いることで，GPU 上した画像を検索する．その処理では，まずクエリ画像から特徴

で特徴ベクトル毎のスレッド並列で処理が可能である．atomi-

量 Q を抽出し，GPU 上のグローバルメモリに転送する．そ

cAdd とは，排他制御を持つ加算関数であり，同一アドレスへして，GPU 上で LSH が適用され，ハッシュ値に変換される．

の書き込みの衝突を回避できる．これを用いて出現頻度集計用

LSH による圧縮後，ハッシュ値とデータベースとの照合を行の配列要素をインクリメントすることにより，最頻出の画像 ID

う．最後に，照合の結果が検索結果となる．

を検出する．しかし，atomicAdd による書き込みは異なるス

クエリ処理のうち，GPU 上で行われる処理は，特徴ベクト

レッドによる同一アドレスへのアクセスが制限されるため，低ルの圧縮とデータベースとの照合の二つに大別される．以下で速な処理の原因となる．シェアードメモリ上のアドレスに対しそれぞれの処理について説明する．

て atomicAdd の書き込みを行う事により，比較的高速な処理

4. 2. 1 特徴ベクトルの圧縮が可能であるが，扱う画像の数が多い場合，メモリ容量が不足

特徴ベクトルの圧縮では，GPU 上のグローバルメモリに転する．

送されたクエリの各特徴ベクトルをハッシュ値へ変換する．一

そこで，本研究では複数の配列を用いて照合処理を高速に行

つのクエリ画像からは数千程度の 128 次元特徴ベクトル Q がう方法を提案する．本研究で提案する照合処理は，一致リスト

抽出されるが，それぞれのベクトルの LSH による変換は独立の取得と最頻出 ID の検出の二つの処理に分けられる．一致リに処理が可能である．また，一つの特徴ベクトル qi ∈ Q の変ストの取得は，クエリから得られたハッシュ値に一致するハッ

換に着目すると，2. 1. 3 節で述べたハッシュ値 hj (qi ) の計算もシュ値をデータベースから検出し，それに対応する全ての画像

それぞれ独立に処理可能である．そこで，各特徴ベクトル qi

ID のリスト (以下，一致リスト) を得る処理である．また，最の変換毎にブロック並列，さらに各ハッシュ値 hj (qi ) の計算

頻出 ID の検出は，一致リストの中で最頻出の画像 ID (以下，毎にスレッド並列で行う．この変換処理は Algorithm 2 で表さ

最頻出 ID) を検出する処理である．

れる．これにより，GPU の並列処理性能を活かすことができ，

a ) 一致リストの取得

高速な処理が可能である．

一致リストの取得では，Algorithm 3 により一致リスト M を-Algorithm 4
Input: M
Output: id
1: M ′ =sort(M )420 15

2: B[M ′ .size()] = 1
3: for i = 0 : M ′ .size() − 1 do in thread parallel
4:

B[i] = M ′ [i] ! = M ′ [i + 1]

5: C = scan(B)
6: for i = 0 : B.size() − 1 do in thread parallel

図3

配列 P，N，F を用いた配列 M の作成

7:

if B[i] == 1 then

8:

O[C [i]] = i

9: S [0] = O[0] + 1
10: for i = 1 : O.size() − 1 do in thread parallel11:

S [i] = O[i] − O[i − 1]

12: r = reduction(S )13: id = M ′ [O[r ]]M の中で最頻出の画像 ID を検出する．その検出は，M をソートした配列 M ′ を用いて行われる．この時，M ′ の値が変化す
るインデックスを境界とし，境界間の部分配列を区間とした時，
図 4 配列 M に対する，配列 M’，B，C，O，S の作成取得する．一致リストは適切な id val の値を配列 M にコピーすることによって取得可能である．GPU のブロック並列で Mを取得するために，各ブロックがコピーする id val の値の先頭
インデックスとコピーサイズ及び書き込み先の先頭インデック
スを事前に計算しておく．以降は，これらを格納した配列をそれぞれ P , N , F とする．さらにこれらの配列を用いて M を効率的に作成する．以下で，配列 P , N , F の作成と，P , N , F を用いた M の作成の詳細を説明する．

P , N , F の作成では，2 行目の処理で各ハッシュ値に対して
f val の中から一致するハッシュ値を二分探索で検索し，マッチした配列のインデックスを index とする．このとき，二分探索がマッチしなかった場合，index は −1 とする．次に，3〜6行目の処理と 7〜10 行目の処理で，配列 P と N を作成する．
これらの処理は，各クエリのハッシュ値に対してスレッド並列で処理を行う．最後に，11 行目の処理で N に対して scan プリ
ミティブを適用し，配列 F を作成する．これらの配列はそれぞれ，P はハッシュ値にマッチする id ptr の値，N はハッシュ値にマッチする ID の個数，F は F [i] が i 番目以前のハッシュ値にマッチする ID の個数を表す．また，図 2 は，これらの配列の作成例を示している．
配列 M は，配列 P , N , F をもとに，12 行目の処理によって作成される．これは，id val の P [i] 番目から N [i] 個の値を Mの F [i] 番目から N [i] 個の要素としてコピーすることを表している．このコピーを各 i に対してブロック並列，コピーする各値に対してスレッド並列で行う．この時，配列 F によって書き込み先の先頭インデックスとサイズが事前に判明しているため書き込みの衝突が起こらず，高速な並列処理が可能である．図

3 は，このコピーの例である．
b ) 最頻出 ID の検出
最頻出 ID の検出では，一致リストの取得によって得られたサイズが最大となる区間の要素が最頻出の画像 ID である．本
手法では，Algorithm 4 によりこの検出処理を高速に行う．
まず，1 行目では，M をソートし，M ′ とする．次に，2〜4行目の処理により，M ′ の境界位置を表す配列 B を作成する．
その後，5 行目の処理で，配列 B に scan プリミティブを適用
し，配列 C を作成する．これは，M ′ のそれぞれの位置の区間
番号を表している．さらに，6〜9 行目の処理で，各区間の終点までのオフセットを表す配列 O を作成する．最後に，10〜11
行目の処理で，各区間のサイズを表す配列 S を作成する．これ
らの処理は，それぞれの配列作成の各要素の計算が独立に行えるため，スレッド並列での処理が可能である．加えて，sort やscan といった GPU 向けのライブラリが存在するプリミティブを利用している．そのため，GPU 上での高速な処理が可能である．図 4 は，これらの配列の作成例である．
ここで，サイズが最も大きくなる区間に属している画像 ID
が検索結果となる．そこで，サイズが最大となる区間番号を検出するために，reduction プリミティブを適用する．しかし，

GPU 向けのライブラリには値が最大値をもつような配列インデックスを検出する関数は存在しないため，単純な reduction
処理を拡張した関数を実装した．この関数を利用して，12 行
目の処理により区間のサイズが最大となる区間番号 r を検出する．このとき，13 行目で表される処理により，検索結果が求め
られる．

5. 実

験

本節では，評価実験について説明する．提案手法は既存手法と比較して高速な検索を実現するが，LSH を用いてデータを圧縮しているため，精度が低下することが予想される．そこで，提案手法の検索速度に加えて検索精度を検証する．

5. 1 データセット
本実験では，公開データセットと独自データセットを含む三つのデータセットを用いた．公開データセットとして，ukbench

と First-Person Social Interactions Dataset を用いた．また，

18,000

独自データセットは，ビデオカメラを用いて作成した．以下に

16,000

それぞれ詳細を示す．

14,000

セットであり，2,550 個の被写体をそれぞれ異なる角度から 4枚ずつ撮影された写真で構成されている．それぞれの画像の解像度は 640×480 であり，約 15 百万本の SIFT 特徴ベクトルが抽出された．

LSH-CPU
LSH-GPU

12,000
検索速度⽐

ukbench [13] （注 4） は，10,200 枚の画像からなる画像データ

10,000
8,000
5,329

6,000
4,000
2,000

（注 5）

First-Person Social Interactions Dataset [10](以下，FP-

1,095

0

2,070

1,209

1

1

ukbenchSID) は，114 個の動画からなるデータセットであり，この動画の一つである AlinDay1-001 の各フレームを画像集合として用

16,562

SIFT-CPU

783

1

FPSID

UTWD

図 5 SIFT-CPU に対する各検索手法の検索速度比．

いた．これは，フレーム数 31,277 枚，解像度 1280×720 の動画である．これを，解像度 640×360 にリサイズして用いて，約43 百万本の SIFT 特徴ベクトルが抽出された．
ビデオカメラを装着し筑波大学構内を歩きながら撮影した
位置情報付きの独自データセット (以下，UTWD) も利用した．UTWD は OLYMPUS STYLUS TG-Tracker を胸部前方に固定した状態で筑波大学構内を歩くことにより撮影した．本
実験では，平砂学生宿舎共用棟前 (36.097688, 140.103352 付近) から 2C 棟 (36.111261, 140.100812 付近) までの経路 (以下，UTWD-long) と，大学会館前 (36.105397, 140.102480 付近) から中央図書館前 (36.109306, 140.101669 付近) まで経路

(以下，UTWD-short) を撮影し，これらの各フレームを画像集合として利用した．なお，UTWD-short の経路は UTWD-longの経路中に含まれる．UTWD-long はフレーム数 81,936 枚，解像度 1280×720 の動画である．これを，解像度 640×360 にリサイズして用いて，約 96 百万本の SIFT 特徴ベクトルが抽出された．一方，UTWD-short はフレーム数 25,980 枚，解像度

1280×720 の動画である．また，UTWD では歩行者の位置推定に利用するために動画の各フレームに TG-Tracker で得られる GPS 情報を付加した．ただし，TG-Tracker では 2 秒毎のGPS 情報のみが得られるため，直接位置情報を得ることができない時刻については線形補間を用いておおよその位置情報を求めた．なお，このデータセットは一般の歩行者の映像が含まれているため非公開である．

5. 2 比 較 手 法
実験の比較手法として，提案手法を含めて三つの画像検索プ
ログラムを実装した．三つの手法は以下の通りである．

SIFT-CPU: 2. 1. 2 節で述べた検索手法の並列実行プログ
ラム．

LSH-CPU: 2. 1. 3 節で述べた検索手法の並列実行プログ
ラム．

LSH-GPU: 4. 節で述べた提案手法のプログラム．
プログラムは C++及び CUDA で記述し，nvcc 7.0 を用いてコンパイルを行った．また，最適化オプションとして -O3 を用いた．CPU での並列化には OpenMP（注 6） を用い，SIFT 特（注 4）：Recognition Benchmark Images www.vis.uky.edu/~stewe/ukbench/（注 5）：First-Person Social Interactions Dataset http://ai.stanford.

edu/~alireza/Disney/（注 6）：OpenMP www.openmp.org/

徴ベクトルの抽出には OpenCV（注 7） を用いた．
プログラムは，CentOS release 6.7 を搭載しているマシン上にて実行した．このマシンは，Intel Xeon Processor E5-2687W

v2 と，NVIDIA Tesla K40 を搭載している．メインメモリは128 GB，デバイスメモリは 12 GB である．
5. 3 実 験 方 法
前述のデータセットとプログラムを用いて，以下の二つの実験を行った．

a) 実

験1

各プログラムにおける検索時間を比較するために，ukbenchと FPSID 及び UTWD-long を画像データベースとして，複数のクエリ画像に対して各プログラムで検索を行い平均検索時間を計測した．ukbench と FPSID では画像データベース中に含まれる画像をクエリ画像とし，UTWD-long に対しては

UTWD-short に含まれる画像各 100 件をクエリ画像とした．ただし，クエリ処理の内，特徴ベクトルを抽出してから画像 IDを返すまでの処理を検索時間として計測した．

b) 実

験

2

LSH による精度変化を検証するために，非圧縮の SIFT 特
徴ベクトルを用いたプログラムと LSH を用いたプログラムにおける検索結果の位置情報に基づく正解率を比較した．このとき，LSH を用いたプログラムについてはパラメータ W とハッ
シュ値サイズ k をそれぞれ変化させて正解率を計算した．画像
データベースとして UTWD-long を用い，クエリ画像として

UTWD-short 中の画像 100 件を用いた．
検索結果はクエリ画像が撮影された地点と検索結果の画像の地点が 15m 以内であれば正解とした．これは位置情報として

GPS を用いていることによる誤差を考慮し，実際のアプリケー
ションにおける許容誤差より大きく設定している．このとき，

GPS 情報から距離への変換には Hubeny の公式 [18] を用いた．
Hubeny の公式は緯度経度情報から距離を計算可能な式であり，短距離であれば高精度な距離を取得可能である．

5. 4 実 験 結 果
本節では，実験 1 と実験 2 の結果および考察を順に示す．

a) 実

験

1

まず，実験 1 の結果を図 5 に示す．図 5 では，検索時間の逆（注 7）：OpenCV | OpenCV opencv.org

1

1
0.8

0.8
0.62
0.49

0.62

0.56

0.6
正解率

正解率

0.6

0.59

0.4

0

0.2

0.15

0.2

0.06

0
10

50

0
100

0.4

0.4

200

300

400

500

0

0.12
0
8

1000

16

64

0

128

256

ハッシュ値サイズ k

パラメータ W

図 6 異なる W における正解率．

数を検索速度とし，各データセットにおける SIFT-CPU の検

32

0.02

図7

異なる k における正解率．

可能であることを示した．

索速度を 1 としたときの各検索手法の検索速度比を示してい

今後の展望としては，クエリ処理のさらなる高速化のためにる．横軸がデータセットと検索手法であり，縦軸に検索速度比

特徴ベクトルの高速な抽出方法を検討する必要がある．また，を示している．ただし，ハッシュ値サイズは k = 16，パラメー

精度向上のために，LSH における適切なパラメータやハッシュタは W = 300 に固定とした．UTWD データセットの結果に

値サイズの検討も必要である．

着目すると，SIFT-CPU と比較して約 1.7 万倍，LSH-CPU と比較して約 3.1 倍と高速である．一方，ukbench データセットと FPSID データセットの結果に着目すると，UTWD データ

謝辞

本研究は JSPS 科研費 JP26280037 の助成を受けたものです．

セットよりも速度が向上していないことが分かる．この原因として，ukbench と FPSID のクエリ画像あたりの特徴ベクトル
数が少ないことが考えられる．クエリ画像あたりの特徴ベクト
ル数は UTWD では約 3,940 に対して，ukbench と FPSID ではそれぞれ約 1,601 と約 1,441 である．そのため，並列度が上がらずに性能向上が小さくなったと考えられる．

b) 実

験

2

図 6 にパラメータ W を変化させた際の正解率を示す．この時，ハッシュ値サイズは k = 16 に設定した．図 6 では，横軸を W とし，縦軸に各 W における正解率を示している．また，
非圧縮 SIFT を用いたときの正解率 (= 0.95) を水平線で表している．これによると，W の値は検索結果の正解率に大きく影響を与えることが分かる．特に，今回の結果の中では，W = 300とした時の正解率が最も高く，62 % の正解率を実現している．
図 7 にハッシュ値サイズ k を変化させた際の正解率を示す．
この時，パラメータは W = 300 に固定して計測した．図 7 では 6 と同様に，横軸を k とし，縦軸に各 k における正解率を示している．W と同様に，k も検索結果の正解率に影響することが分かる．

6. ま と め
本研究では，SIFT 特徴量を用いた画像検索の GPU による高速化手法を提案した．その際，LSH を用いて特徴ベクトルを圧縮することにより，空間コストと時間コストの削減を実現した．また，GPU 上での処理に適したデータ構造およびアルゴ
リズムにより，高速な処理を実現した．
実験により，CPU による検索手法と比較して約 3.1 倍程度の高速化が可能であることを示した．また，実アプリケーショ
ンでの適用を想定した精度評価を行い，62 % の正解率が実現

文

献[1] Jurandy Almeida, Ricardo da S Torres, and Siome Goldenstein. Sift applied to cbir. Revista de Sistemas de Informacao da FSMA n, Vol. 4, pp. 41–48, 2009.
[2] Alexandr Andoni and Piotr Indyk. Near-optimal hashingalgorithms for approximate nearest neighbor in high dimensions. Commun. ACM, Vol. 51, No. 1, pp. 117–122, 2008.
[3] Adrien Auclair, Laurent D Cohen, and Nicole Vincent. How
to use sift vectors to analyze an image with database templates. In Adaptive Multimedia Retrieval: Retrieval, User,and Semantics, pp. 224–236. Springer, 2007.
[4] R. Baeza-Yates and B. Ribeiro-Neto. Modern Information
Retrieval: The Concepts and Technology Behind Search,
Vol. 2. Addison Wesley, 2011.
[5] Lee Carraher, Philip A. Wilsey, and Fred S. Annexstein.
A GPGPU algorithm for c-approximate r-nearest neighborsearch in high dimensions. In 2013 IEEE International
Symposium on Parallel & Distributed Processing, Workshops and Phd Forum, Cambridge, MA, USA, May 20-24,
2013, pp. 2079–2088, 2013.
[6] Ali Cevahir and Junji Torii. Gpu-enabled high performance
online visual search with high accuracy. In 2012 IEEE International Symposium on Multimedia, ISM 2012, Irvine,
CA, USA, December 10-12, 2012, pp. 413–420, 2012.
[7] Vijay Chandrasekhar, Mina Makar, Gabriel Takacs, David
Chen, Sam S Tsai, Ngai-Man Cheung, Radek Grzeszczuk,
Yuriy Reznik, and Bernd Girod. Survey of sift compression schemes. In Proc. Int. Workshop Mobile Multimedia
Processing, pp. 35–40. Citeseer, 2010.
[8] Jian Cheng, Cong Leng, Jiaxiang Wu, Hainan Cui, andHanqing Lu. Fast and accurate image matching with cascade hashing for 3d reconstruction. In 2014 IEEE Conference on Computer Vision and Pattern Recognition, CVPR
2014, Columbus, OH, USA, June 23-28, 2014, pp. 1–8,
2014.
[9] Mayur Datar, Nicole Immorlica, Piotr Indyk, and Vahab S.
Mirrokni. Locality-sensitive hashing scheme based on pstable distributions. In Proceedings of the 20th ACM Sym-

[10]

[11]

[12]
[13]

[14]

[15]

[16]

[17]

[18]

posium on Computational Geometry, Brooklyn, New York,USA, June 8-11, 2004, pp. 253–262, 2004.
Alireza Fathi, Jessica K. Hodgins, and James M. Rehg. Social interactions: A first-person perspective. In 2012 IEEE
Conference on Computer Vision and Pattern Recognition,
Providence, RI, USA, June 16-21, 2012, pp. 1226–1233,
2012.
Piotr Indyk and Rajeev Motwani. Approximate nearest
neighbors: Towards removing the curse of dimensionality.
In Proceedings of the Thirtieth Annual ACM Symposiumon the Theory of Computing, Dallas, Texas, USA, May
23-26, 1998, pp. 604–613, 1998.
David G. Lowe.
Object recognition from local scaleinvariant features. In ICCV, pp. 1150–1157, 1999.
David Nistér and Henrik Stewénius. Scalable recognition
with a vocabulary tree. In 2006 IEEE Computer Society
Conference on Computer Vision and Pattern Recognition(CVPR 2006), 17-22 June 2006, New York, NY, USA, pp.
2161–2168, 2006.
John D Owens, Mike Houston, David Luebke, Simon Green,
John E Stone, and James C Phillips. Gpu computing. Proceedings of the IEEE, Vol. 96, No. 5, pp. 879–899, 2008.
Josef Sivic and Andrew Zisserman. Video google: A
text retrieval approach to object matching in videos. In
9th IEEE International Conference on Computer Vision(ICCV 2003), 14-17 October 2003, Nice, France, pp. 1470–
1477, 2003.
W. Weihong and W. Song. A scalable content-based image
retrieval scheme using locality-sensitive hashing. In 2009 International Conference on Computational Intelligence andNatural Computing, Vol. 1, pp. 151–154, June 2009.
亀田能成, 大田友一. 街中での歩行者カメラによるオンライン
位置推定のための検討. 画像の認識・理解シンポジウム MIRU,
pp. 364–369, 2010.
日本測量協会創立 30 周年記念「現代測量学」出版委員会. 現代
測量学, 第 4 巻. 日本測量協会, 1988.