コンピュータビジョンセミナー vol.4~Visual SLAM / SfMで行われる局所特徴量計算のCUDA高速化~(2024/3/27)

7.3K Views

March 28, 24

スライド概要

自動移動ロボットに欠かせない Visual SLAM や、ドローンなどの撮像データから3Dモデルを構築する SfM では、局所特徴量計算が行われています。
当社がGithub で公開した cuda-efficient-features には、局所特徴量計算アルゴリズムである BAD と HashSIFT の実装が含まれます。

本セミナーでは、HashSIFTのアルゴリズムとCUDA高速化について解説していきます。

当社ブログ記事(以下URL)としても公開していますので、是非ご覧ください。

<出典:当社ブログ記事>
・CUDAによる局所特徴量計算の高速化とソースコード公開:
  https://proc-cpuinfo.fixstars.com/2023/12/cuda-efficient-features/
・特徴記述アルゴリズムHashSIFTのCUDA高速化:
  https://proc-cpuinfo.fixstars.com/2023/12/cuda-hash-sift/


<コンピュータビジョンシリーズの過去資料>
・vol.1 OpenCV活用(OpenCVでCUDAを活用するためのGpuMat解説):
  https://www.docswell.com/s/fixstars/ZRXQ72-20220805
・vol.2 視差計算ライブラリ libSGM のアルゴリズム解説と CUDA高速化:
  https://www.docswell.com/s/fixstars/5ENE7J-20221220
・vol.3 視差計算を利用した障害物検出:
  https://www.docswell.com/s/fixstars/ZQ8447-20231027
・vol.4 Visual SLAM / SfMで行われる局所特徴量計算のCUDA高速化
  https://www.docswell.com/s/fixstars/ZM1DJ8-20240327

・フィックスターズの画像処理アルゴリズム開発支援:
  https://www.fixstars.com/ja/services/image-processing

profile-image

フィックスターズは、コンピュータの性能を最大限に引き出すソフトウェア開発のスペシャリストです。車載、産業機器、金融、医療など、幅広い分野での開発経験があります。また、ディープラーニングや機械学習などの最先端技術にも力を入れています。 並列化や最適化技術を駆使して、マルチコアCPU、GPU、FPGA、量子アニーリングマシンなど、さまざまなハードウェアでソフトウェアを高速化するサービスを提供しています。さらに、長年の経験から培ったハードウェアの知識と最適化ノウハウを活かし、高精度で高性能なアルゴリズムの開発も行っています。       ・開催セミナー一覧:https://www.fixstars.com/ja/seminar   ・技術ブログ :https://proc-cpuinfo.fixstars.com/

シェア

またはPlayer版

埋め込む »CMSなどでJSが使えない場合

(ダウンロード不可)

関連スライド

各ページのテキスト
1.

コンピュータビジョンセミナーvol.4 Visual SLAM / SfM で行われる局所特徴量計算のCUDA高速化 Copyright © Fixstars Group

2.

本セミナーの構成 1. 会社紹介 2. Visual SLAM / SfM における局所特徴量計算 3. HashSIFTのアルゴリズム 4. HashSIFTのCUDA高速化 Copyright © Fixstars Group 1

3.

本セミナーの位置づけ 弊社でサービス展開しているコンピュータビジョン領域での様々な技術情報を、 コンピュータビジョンセミナーとして発信しています ⚫ Vol.1 OpenCVでCUDAを活用するためのGpuMat解説 ⚫ Vol.2 視差計算ライブラリ libSGM のアルゴリズム解説 ⚫ Vol.3 視差計算を利用した障害物検出 ⚫ Vol.4 Visual SLAM / SfM で行われる 局所特徴量計算の CUDA高速化 ⚫ こんな方に向いています ○ ○ 局所特徴量計算のアルゴリズムを知りたい 実用的なアルゴリズムの GPU 高速化事例を知りたい Copyright © Fixstars Group 2

4.

発表者紹介 冨田 明彦 高木 章洋 ソリューションカンパニー 営業企画 ソリューション第二事業部 エグゼクティブエンジニア 2008年に入社。金融、医療業界において、ソ フトウェア高速化業務に携わる。その後、新規 事業企画、半導体業界の事業を担当し、現職。 2014年に新卒入社。自動運転向けの画像認識 アルゴリズムの開発、高速化などに携わる。 OSS活動として、libSGMの開発にも携わる。 Copyright © Fixstars Group 3

5.

フィックスターズの ご紹介 Copyright © Fixstars Group 4

6.

フィックスターズの強み コンピュータの性能を最大限に引き出す、ソフトウェア高速化のエキスパート集団 ハードウェアの知見 アルゴリズム実装力 各産業・研究分野の知見 目的の製品に最適なハードウェアを見抜き、 その性能をフル活用するソフトウェアを開 発します。 ハードウェアの特徴と製品要求仕様に合わ せて、アルゴリズムを改良して高速化を実 現します。 開発したい製品に使える技術を見抜き、実 際に動作する実装までトータルにサポート します。 Copyright © Fixstars Group 5

7.

サービス概要 お客様専任のエンジニアが直接ヒアリングを行い、高速化を実現するために乗り越えるべき 課題や問題を明確にしていきます。 高速化のワークフロー お客様 オリジナルソースコードのご提供 高速化したコード コンサルティング 高速化 サポート 先行技術調査 アルゴリズムの改良・開発 レポートやコードへのQ&A 性能評価・ボトルネックの特定 ハードウェアへの最適化 実製品への組込み支援 レポート作成 Copyright © Fixstars Group 6

8.

サービス提供分野 半導体 産業機器 金融 自動車 ● NAND型フラッシュメモリ向け ファームウェア開発 ● 次世代AIチップの開発環境基盤 生命科学 ● Smart Factory実現への支援 ● マシンビジョンシステムの高速化 ● 自動運転の高性能化、実用化 ● ゲノム解析の高速化 ● 次世代パーソナルモビリティの 研究開発 ● 医用画像処理の高速化 Copyright © Fixstars Group ● デリバティブシステムの高速化 ● HFT(アルゴリズムトレード)の高速化 ● AI画像診断システムの研究開発 7

9.

サービス領域 様々な領域でソフトウェア高速化サービスを提供しています。大量データの高速処理は、 お客様の製品競争力の源泉となっています。 組込み高速化 GPU向け高速化 AI・深層学習 画像処理・ アルゴリズム開発 FPGAを活用した システム開発 分散並列システム開発 量子コンピューティング 自動車向け フラッシュメモリ向け ソフトウェア開発 ファームウェア開発 Copyright © Fixstars Group 8

10.

画像処理アルゴリズム開発 高速な画像処理需要に対して、経験豊富なエンジニアが 責任を持って製品開発をご支援します。 お客様の課題 ご支援内容 高度な画像処理や深層学習等のアルゴリズム を開発できる人材が社内に限られている アルゴリズム調査・改変 課題に合ったアルゴリズム・実装手法を調査 製品実装に向けて適切な改変を実施 機能要件は満たせそうだが、ターゲット機器 上で性能要件までクリアできるか不安 深層学習ネットワーク精度の改善 様々な手法を駆使して深層学習ネットワークの精度を改善 製品化に結びつくような研究ができていない 論文調査・改善活動 論文調査から最先端の手法の探索 性能向上に向けた改善活動を継続 Copyright © Fixstars Group 9

11.

GPU向け高速化 高性能なGPUの本来の性能を十分に引き出し、 ソフトウェアの高速化を実現します。 お客様の課題 ご支援内容 GPUで計算してみたが期待した性能が出ない GPU高速化に関するコンサルティング GPU/CPUを組み合わせた全体として最適な設 CPU・GPU混在環境でのシステム設計 計がしたい アルゴリズムのGPU向け移植 原価を維持したまま機能を追加するため、も う少し処理を速くしたい GPUプログラム高速化 品質確保のため、精度を上げたく演算量は増 えるが性能は維持したい Copyright © Fixstars Group 継続的な精度向上 10

12.

1. Visual SLAM / SfM における局所特徴量計算 Copyright © Fixstars Group 11

13.

局所特徴量計算 ⚫ 画像の局所的な特徴を検出・記述するための計算 ⚫ Visual SLAMやSfMでは、特に画像間の特徴点マッチングで使用 ⚫ 特徴点検出と特徴量記述から構成される ○ 特徴点検出:画像からマッチングに適した点を検出 (典型的にはコーナー) ○ 特徴量記述:特徴点周辺の局所領域から、回転・スケール変動・照明変動などの 様々な見えの変化に頑健な、特徴ベクトルを抽出 特徴点マッチングの例 ([1]より引用) Copyright © Fixstars Group 12

14.

代表的な局所特徴量計算手法 ⚫ SIFT (Scale-Invariant Feature Transform) [2] ○ DoG画像によるスケールスペースを用いたスケール不変な特徴点の検出 ○ 特徴点の周辺領域から勾配方向ヒストグラムを作成し、特徴量化 ⚫ ORB (Oriented FAST and Rotated BRIEF) [1] ○ 画像ピラミッドを用いたマルチスケールのFAST特徴点検出 ○ 特徴点の周辺領域から選んだ2点の輝度値の大小に1/0を割り当て、バイナリ特徴量化 ○ 省メモリかつ類似度計算を高速に実行できるメリット SIFT特徴量の計算 ORBのサンプリングパターン ([1]より引用) Copyright © Fixstars Group 13

15.

BADとHashSIFTの登場 [3] (Suárez et al. 2021) ⚫ 前述の局所特徴量をベースに、深層学習の学習テクニックを導入し性能向上 ○ Triplet Ranking Loss, Hard Negative Mining, anchor swap, Data Augmentation ⚫ BAD (Box Average Difference) ○ ORBと同じく、輝度差を利用した局所特徴量 ○ 2つの正方領域(Box)の平均輝度値の差を閾値処理してバイナリ特徴量化 ⚫ HashSIFT ○ SIFTと同じく、勾配方向を利用した局所特徴量 ○ SIFT特徴量に線形変換と2値化を適用してバイナリ特徴量化 SIFT特徴量 𝐟 𝐱 𝐁 sgn b11 ⋮ bD1 BADの計算のイメージ ([4]より引用) ⋯ ⋱ ⋯ b1K+1 ⋮ bDK+1 f1 𝐱 f2 𝐱 f3 𝐱 ⋮ fK 𝐱 1 𝐡 𝐱 +1 −1 → +1 ⋮ ⋮ −1 𝐃 𝐱 d 1 𝐱 1 d2 𝐱 0 → 1 = d3 𝐱 ⋮ ⋮ ⋮ ⋮ 0 dD 𝐱 HashSIFTの計算のイメージ Copyright © Fixstars Group 14

16.

BADとHashSIFTのマッチング性能 ⚫ ホモグラフィ推定によるマッチング性能比較※ ○ 2枚の画像から特徴点を検出し、評価対象の特徴量を抽出 ○ 特徴点マッチングを行い、その結果からホモグラフィ行列を推定 (RANSAC使用) ○ RANSACで得られたインライア数を比較 ⚫ BADとHashSIFTのマッチング性能が高いことが確認できる ORB インライア数:35 BAD インライア数:128 Copyright © Fixstars Group SIFT インライア数:64 HashSIFT インライア数:183 ※論文著者が公開しているサンプルコード[5]を元に作成 15

17.

BADとHashSIFTの処理時間 ⚫ モバイルプロセッサ上でも高速 ⚫ 優れた精度と計算コストのバランス ORB系 SIFT系 CNN系 特徴量計算の処理時間比較 ([3]より引用) Oxfordデータセット[11]の画像から検出した特徴点 (1画像最大2000点)に対する、特徴量記述の平均処理時間 特徴量計算の精度と計算コストのトレードオフ ([3]より引用) ※HashSIFTがSIFT(OpenCV実装)より速い理由 OpenCV実装では事前にガウシアンピラミッドを作成し、スケールレベル毎に処理するのに対し、HashSIFTは入力画像を直接処理するだけのため Copyright © Fixstars Group 16

18.

モチベーション ⚫ 高速・軽量な局所特徴量計算はSLAM・SfMにおいて有用 ○ SLAMではリアルタイム性が求められる ○ SfMでは大量の画像を処理する ⚫ BAD・HashSIFTはこれらの期待を満たす有力手法だが、 GPUでさらに高速化したい ○ 特徴量記述は特徴点毎に並列に計算できるため、高速化の期待度大 ⚫ 本セミナーではHashSIFTにフォーカスして解説 Copyright © Fixstars Group 17

19.

2. HashSIFTのアルゴリズム Copyright © Fixstars Group 18

20.

HashSIFT計算の流れ ⚫ 局所画像パッチ作成 ○ 画像中の特徴点周りの領域を、特徴点のスケールとオリエンテーション※1によって切り出し、 アフィン変換によって32×32ピクセルのパッチに変換 ⚫ SIFT特徴量計算 ○ パッチから勾配方向ヒストグラムを作成し、128次元のSIFT特徴量を計算 ⚫ 線形変換+2値化 ○ SIFT特徴量を線形変換でD次元※2のベクトルに変換後、閾値処理によりバイナリ特徴量化 sgn 局所画像パッチ 輝度勾配 勾配方向ヒストグラム SIFT特徴量 b11 ⋮ bD1 ⋯ ⋱ ⋯ b1K+1 ⋮ bDK+1 f1 𝐱 f2 𝐱 f3 𝐱 ⋮ fK 𝐱 1 d1 𝐱 +1 1 d2 𝐱 −1 0 → +1 → 1 = d3 𝐱 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ 0 −1 dD 𝐱 線形変換+2値化 ※1 特徴点検出およびスケール・オリエンテーションの計算は、SIFTベースやORBベースなど、何を採用しても良い ※2 D = 256 or 512 Copyright © Fixstars Group 19

21.

SIFT特徴量計算 (1/2) ⚫ パッチ𝐼の微分画像𝐼𝑥 , 𝐼𝑦 を作成し、勾配強度𝑚と勾配方向𝜃を計算 ○ 𝐼𝑥 𝑥, 𝑦 = 𝐼 𝑥 + 1, 𝑦 − 𝐼 𝑥 − 1, 𝑦 ൝ 𝐼𝑦 𝑥, 𝑦 = 𝐼 𝑥, 𝑦 + 1 − 𝐼 𝑥, 𝑦 − 1 ○ 𝑚 𝑥, 𝑦 = ○ 𝜃 𝑥, 𝑦 = tan−1 𝐼𝑥 𝑥, 𝑦 2 + 𝐼𝑦 𝑥, 𝑦 2 𝐼𝑦 𝑥,𝑦 𝐼𝑥 𝑥,𝑦 局所画像パッチ (32×32) Copyright © Fixstars Group 輝度勾配 (θにより色分け) 20

22.

SIFT特徴量計算 (2/2) ⚫ 局所ヒストグラム計算 ○ パッチを8×8ピクセルの小ブロックで分割 (ブロック数は4×4) ○ ブロック毎に8方向の勾配方向ヒストグラムを作成 ■ 対応するビンに、パッチ中心からの距離で重みづけした勾配強度𝑚を累積 ⚫ 特徴ベクトル化+正規化 ○ ヒストグラムを1次元化し、4×4×8=128次元の特徴ベクトル化 ○ ノルムが1となるよう特徴ベクトルを正規化 輝度勾配 勾配方向ヒストグラム Copyright © Fixstars Group SIFT特徴量 (128次元) 21

23.

線形変換+2値化 ⚫ SIFT特徴量を線形変換でD次元のベクトルに変換 ○ 変換行列の係数は事前に学習(後述) ○ n個の特徴点全てのSIFT特徴量をまず計算し、1つの行列に結合後、行列-行列積を行う ⚫ 各要素を閾値処理によって、1/0のバイナリ特徴量に変換 φ ステップ関数 φ x =ቊ 0 x≤0 1 x>0 b11 ⋮ bD1 ⋯ b1K+1 ⋱ ⋮ ⋯ bDK+1 変換行列 f1 𝐱1 f2 𝐱 1 f3 𝐱 1 ⋮ fK 𝐱 1 1 f1 𝐱2 f2 𝐱 2 f3 𝐱 2 ⋮ fK 𝐱 2 1 f1 𝐱n f2 𝐱 n ⋯ f3 𝐱 n ⋮ fK 𝐱 n 1 SIFT特徴量 (結合後) 末尾行の1はバイアス項 Copyright © Fixstars Group d1 𝐱1 d2 𝐱1 = d3 𝐱1 ⋮ ⋮ dD 𝐱1 d1 𝐱2 d2 𝐱2 d3 𝐱2 ⋮ ⋮ dD 𝐱2 d1 𝐱n d2 𝐱n ⋯ d3 𝐱n ⋮ ⋮ dD 𝐱n HashSIFT特徴量 (結合後) 22

24.

HashSIFTの学習 (1/3) ⚫ Liberty dataset[9]を用いた学習 ○ ラベル付けされたパッチ画像の集合 ⚫ パッチのtriplet(三つ組)の集合𝒯を作成 ○ 𝒯= 𝐚i , 𝐩i , 𝐧i N i=1 ⚫ Triplet Ranking Loss(TRL)による目的関数𝓛𝐁 T 学習用パッチの例 ([9]より引用) はSIFT特徴量を変換した、+1または-1から構成されたベクトル ○ 𝐃 𝐱 = sgn 𝐁 𝐟 𝐱 T , 1 ○ ෩ 𝐱 = tanh 𝐁 𝐟 𝐱 T , 1 sgnをtanhで近似し、𝐃 ○ T෩ T෩ ෩ ෩ 𝓛𝐁 = σN i=1 τ − 𝐃 𝐚i 𝐃 𝐩i + 𝐃 𝐚i 𝐃 𝐧i 𝐚i と𝐩i の類似度が大きくなるとLossが減少 T とすると + 𝐚i と𝐧i の類似度が小さくなるとLossが減少 学習パッチのtriplet 𝐚i :anchorパッチ 𝐩i :positiveパッチ (𝐚i と同一ラベルのパッチ.正例) 𝐧i :negativeパッチ (𝐚i と異なるラベルのパッチ.負例) Copyright © Fixstars Group tanh関数 23

25.

HashSIFTの学習 (2/3) ⚫ 確率的勾配降下法による𝓛𝐁 の最小化 𝐚i , 𝐩i , 𝐧i N i=1 の部分集合𝒫をサンプリングして勾配を計算 ○ 𝒯= ○ 𝐁の要素を繰り返し更新 ⚫ Hard Negative Mining(HNM)[7]とanchor swap[6]によるサンプリング ○ HNM:現時点で間違えやすい(=パッチ間距離が近い)負例を𝒫に含める ○ anchor swap:dist 𝐩i , 𝐧i < dist 𝐚i , 𝐧i ならば、 𝐚i と𝐩i を入れ替える (よりHardに) ⚫ Data Augmentation ○ パッチに様々な変形やノイズを付加 ○ 過学習を抑制 Hard Negative Mining ([7]より引用) Copyright © Fixstars Group 24

26.

HashSIFTの学習 (3/3) ⚫ Binary Classification Lossによる学習と比較して、 TRL,HNM, Data Augmentationのそれぞれに改善効果を確認 HPatches Image Matching Task[11]における、HashSIFTのAblation Study ([3]より引用) Copyright © Fixstars Group 25

27.

3. HashSIFTのCUDA高速化 Copyright © Fixstars Group 26

28.

CUDAの簡単なおさらい ⚫ スレッドおよびメモリの階層 スレッドの階層構造 • スレッド間に階層構造がある • 近いスレッド同士はより密に通信・同期を行うことができる メモリの階層構造 • おおむねスレッドの階層構造と対応 Fixstars CUDA高速化セミナーvol.1 ~画像処理アルゴリズムの高速化~ より引用 https://www.docswell.com/s/fixstars/K24MYM-20220527 Copyright © Fixstars Group 27

29.

HashSIFTのCUDA化の方針 ① 特徴点数分のSIFT計算を行うCUDAカーネルの実装 ② 特徴点数分の行列積はcuBLASライブラリを利用 ③ 特徴点数分の2値化を行うCUDAカーネルの実装 φ ①SIFT特徴量計算×n ③2値化 Copyright © Fixstars Group b11 ⋮ bD1 ⋯ b1K+1 ⋱ ⋮ ⋯ bDK+1 f1 𝐱1 f2 𝐱1 f3 𝐱1 ⋮ fK 𝐱1 1 f1 𝐱2 f2 𝐱 2 f3 𝐱 2 ⋮ fK 𝐱 2 1 ⋯ f1 𝐱n f2 𝐱 n f3 𝐱 n ⋮ fK 𝐱 n 1 =D ②行列積 28

30.

SIFT特徴量計算のCUDA実装 Copyright © Fixstars Group 29

31.

スレッドの割り当て方の検討 ① 1特徴点に1スレッドを割り当て ○ 1特徴点内のSIFT計算はシーケンシャルに実行 ○ スレッド間の同期や競合を気にする必要がなく、実装が簡単 ② 1特徴点に1ブロックを割り当て ← 今回採用 ○ 1特徴点内のSIFT計算に複数スレッドを割り当てることができる ○ Shared Memoryによる高速なデータの共有が可能 ○ スレッド間の同期や競合をケアする必要がある 1特徴点に1ブロックを割り当て Copyright © Fixstars Group 30

32.

1ブロック当たり何スレッドにするか? ① 1ブロック-32×32スレッド ○ パッチに対する処理を最大限並列化できるが、正規化処理でスレッドが余る ○ スレッド間の同期コスト大 ② 1ブロック-32スレッド (1Warp) ○ スレッドの余りはないが、並列度は低下 ○ スレッド間の同期コスト低 ③ 1ブロック-32×4スレッド ← 今回採用 ○ ①と②のバランスをとった選択 32×32要素への処理 Copyright © Fixstars Group 128要素への処理 31

33.
[beta]
実装例:カーネル関数およびカーネル関数呼び出し
__global__ void computeSIFTKernel(const PtrStepSzb image,
const PtrStepSz<KeyPoint> keypoints, PtrStepf responses)
{
// パッチ,ヒストグラム,特徴ベクトルをShared Memoryに配置
__shared__ uchar patch[PATCH_H][PATCH_W];
__shared__ float histogram[R_BINS + 2][C_BINS + 2][ORI_BINS + 2];
__shared__ float descriptor[DESCRIPTOR_SIZE];
// 特徴点のID (0 ~ 特徴点数-1)
const int i = blockIdx.x;

void computeSIFTs(const GpuMat& image, const GpuMat& keypoints,
GpuMat& responses)
{
const dim3 block(32, 4); // 32×4のスレッドを生成
const dim3 grid(keypoints.rows, 1); // 特徴点数のブロックを生成
computeSIFTKernel<<<grid, block>>>(image, keypoints, responses);
CUDA_CHECK(cudaGetLastError());
}

カーネル関数呼び出し

// 局所画像パッチ作成
rectifyPatch(image, keypoints[i], patch);
__syncthreads();
// 勾配方向ヒストグラム作成
createGradientHistogram(patch, histogram);
__syncthreads();

__syncthreads()により、直前の計算結果が
Shared Memory内に全て格納されるまで待つ

// 1次元化 + 正規化
createFeatureVector(histogram, descriptor);
__syncthreads();
// 結果をグローバルメモリに書き込み
for (int k = threadIdx1D(); k < DESCRIPTOR_SIZE; k += blockSize1D())
responses(i, k) = descriptor[k];
}

SIFT計算のカーネル関数

※実際の実装より簡略化しています
Copyright © Fixstars Group

32

34.
[beta]
実装例:局所画像パッチ作成~輝度勾配計算
__device__ void rectifyPatch(const PtrStepSzb image,
const Keypoint& keypoint, uchar patch[PATCH_H][PATCH_W])
{
// キーポイントのスケールとオリエンテーションからアフィン変換行列を作成
const Mat23f M = getAffineTransform(keypoint);
for (int y = threadIdx.y; y < PATCH_H; y += blockDim.y) {
for (int x = threadIdx.x; x < PATCH_W; x += blockDim.x) {

__device__ void createGradientHistogram(uchar patch[PATCH_H][PATCH_W],
float histogram[R_BINS + 2][C_BINS + 2][ORI_BINS + 2])
{
for (int y = threadIdx.y; y < PATCH_H; y += blockDim.y) {
for (int x = threadIdx.x; x < PATCH_W; x += blockDim.x) {
// 中心からの距離を元に投票の重みを計算
const float magScale = expf(distScale * normsq(x - cx, y - cy));

// アフィン変換から入力画像の参照位置を計算
const float u = M(0, 0) * x + M(0, 1) * y + M(0, 2);
const float v = M(1, 0) * x + M(1, 1) * y + M(1, 2);

// 勾配強度と勾配方向の計算
const float dx = patch[y + 0][x + 1] - patch[y + 0][x - 1];
const float dy = patch[y - 1][x + 0] - patch[y + 1][x + 0];
const float mag = magScale * sqrtf(normsq(dx, dy));
const float ori = atan2f(dy, dx);

// 入力画像からサンプリング
patch[y][x] = sampleBilinear(image, u, v);
}

// ヒストグラムへの投票
...

}
}

}

局所画像パッチ作成

}
}

小技:threadIdx初期化+blockDimインクリメント
このようにループを記述することで、
スレッド数 < 要素数な処理に柔軟に対応できる

輝度勾配計算

今回の場合、32×4スレッドで32×32のパッチを、
右図のような順序でアクセスしていく

※実際の実装より簡略化しています
Copyright © Fixstars Group

33

35.

実装例:ヒストグラムへの投票 // 勾配強度と勾配方向の計算 const float mag = magScale * sqrtf(normsq(dx, dy)); const float ori = atan2f(dy, dx); // ブロック座標に変換 const auto [ri, rf] const auto [ci, cf] const auto [oi, of] + = = = 整数部と小数部に分解 separateIF(getRowBin(y)); separateIF(getColBin(x)); separateIF(getOriBin(ori)); 周辺のブロックに投票 c // 小数部に従って投票量(mag)を周辺ブロックに分配 const auto [ v0, v1] = distribute(mag, rf); const auto [ v00, v01] = distribute(v0, cf); const auto [ v10, v11] = distribute(v1, cf); const auto [v000, v001] = distribute(v00, of); const auto [v010, v011] = distribute(v01, of); const auto [v100, v101] = distribute(v10, of); const auto [v110, v111] = distribute(v11, of); // 周辺ブロックに投票 atomicAdd(&histogram[ri atomicAdd(&histogram[ri atomicAdd(&histogram[ri atomicAdd(&histogram[ri atomicAdd(&histogram[ri atomicAdd(&histogram[ri atomicAdd(&histogram[ri atomicAdd(&histogram[ri + + + + + + + + 1][ci 1][ci 1][ci 1][ci 2][ci 2][ci 2][ci 2][ci + + + + + + + + 1][oi 1][oi 2][oi 2][oi 1][oi 1][oi 2][oi 2][oi + + + + + + + + 0], 1], 0], 1], 0], 1], 0], 1], r v000); v001); v010); v011); v100); v101); v110); v111); 複数スレッドが同時にアクセス しないよう、atomicAddを使用 } } } ※実際の実装より簡略化しています Copyright © Fixstars Group 34

36.
[beta]
実装例:特徴ベクトル作成 (1次元化 + 正規化)
__device__ void createFeatureVector(float histogram[…],
float descriptor[DESCRIPTOR_SIZE])
{
// 1次元化
for (int r = threadIdx.y; r < R_BINS; r += blockDim.y)
for (int c = threadIdx.x; c < C_BINS; c += blockDim.x)
for (int k = 0; k < ORI_BINS; k++)
descriptor[(r * R_BINS + c) * ORI_BINS + k] = histogram[r][c][k];
// 正規化
if (threadIdx.y == 0)
{
// 2乗和をWarp内にまとめる
float sumSq = 0.f;
for (int i = threadIdx.x; i < DESCRIPTOR_SIZE; i += WARP_SIZE)
sumSq += squared(descriptor[i]);
// Warp ShuffleでWarp内の総和をとる
for (int mask = 16; mask > 0; mask /= 2)
sumSq += __shfl_xor_sync(0xffffffff, sumSq, mask);
const float scale = 1.f / sqrtf(sumSq);
for (int i = threadIdx.x; i < DESCRIPTOR_SIZE; i += WARP_SIZE)
descriptor[i] = scale * descriptor[i];
}

Warp Shuffleによる総和の実装例 ([13]より引用)

}

特徴ベクトル作成 (1次元化 + 正規化)

※実際の実装より簡略化しています
Copyright © Fixstars Group

35

37.

ヒストグラム作成時のatomicAdd衝突回数削減 ⚫ 複数スレッドが同時に同一ビンにアクセス(衝突)すると待ちが発生 ○ ナイーブ実装では32×4スレッドのまとまりが上~下に移動 ○ スレッド同士で参照位置が近く、周辺ブロックへの投票も含めると衝突頻度が高いはず ⚫ 参照位置が分散するようにアクセスパターンを工夫 ○ 効果は微改善程度 for (int y = threadIdx.y; y < PATCH_H; y += blockDim.y) { for (int x = threadIdx.x; x < PATCH_W; x += blockDim.x) { ナイーブ実装のループの書き方 constexpr int BLOCK_SIZE = 8; constexpr int STRIDE_X = 4; constexpr int STRIDE_Y = CELL_SIZE / 4; ナイーブ実装のアクセスパターン 32スレッドが同時に1ブロックを参照 工夫したアクセスパターン 8スレッドが同時に1ブロックを参照 const int niterations = PATCH_W * PATCH_H / (32 * 4); const int tid = threadIdx1D(); for (int iter = 0; iter < niterations; iter++) { const int y = (tid / BLOCK_SIZE) * STRIDE_Y + (iter / STRIDE_X); const int x = (tid % BLOCK_SIZE) * STRIDE_X + (iter % STRIDE_X); 工夫したループの書き方 色々と試行錯誤してこのようになった Copyright © Fixstars Group 36

38.

SIFT計算以降のCUDA実装 Copyright © Fixstars Group 37

39.
[beta]
行列積
⚫ cuBLASライブラリを使用
void linearProjection(const GpuMat& F, const GpuMat& B, GpuMat& H)
{
const float alpha = 1.f;
const float beta = 0.f;
const cublasOperation_t transa = CUBLAS_OP_T;
const cublasOperation_t transb = CUBLAS_OP_N;
// cuBLASハンドルの初期化
cublasHandle_t handle;
cublasCreate_v2(&handle);
// スカラー値をホスト経由で渡す設定
cublasSetPointerMode_v2(handle, CUBLAS_POINTER_MODE_HOST);
// 行列積の実行 H := alpha*B*F + beta*C
cublasSgemm_v2(handle, transa, transb, B.rows, F.rows, B.cols,
&alpha,
B.ptr<float>(), B.cols,
F.ptr<float>(), F.cols,
&beta,
H.ptr<float>(), H.step.cols);
// cuBLASハンドルの削除
cublasDestroy_v2(handle);
}

Copyright © Fixstars Group

38

40.
[beta]
2値化
⚫ float 8要素を2値化し、1byte(8bit)に変換

__global__ void binarizeDescriptorsKernel(const PtrStepf src, PtrStepSzb dst)
{
const int i = blockDim.y * blockIdx.y + threadIdx.y; // 特徴点のID
const int p = blockDim.x * blockIdx.x + threadIdx.x; // 特徴量のバイト位置

uchar byte = 0;
const float* ptrSrc = src.ptr(i) + p * 8;
byte
byte
byte
byte
byte
byte
byte
byte

|=
|=
|=
|=
|=
|=
|=
|=

(*ptrSrc++
(*ptrSrc++
(*ptrSrc++
(*ptrSrc++
(*ptrSrc++
(*ptrSrc++
(*ptrSrc++
(*ptrSrc++

>
>
>
>
>
>
>
>

0)
0)
0)
0)
0)
0)
0)
0)

<<
<<
<<
<<
<<
<<
<<
<<

7;
6;
5;
4;
3;
2;
1;
0;

dst(i, p) = byte;
}

Copyright © Fixstars Group

39

41.

CUDA APIの呼び出し回数削減 ⚫ CUDA API ○ デバイスメモリの確保・解放やcuBLASハンドルの初期化・削除など ○ カーネル関数が高速になると、相対的にボトルネックになることが多い ⚫ 1度確保したデバイスバッファやハンドルをHashSIFTの実装クラス内に保持、 使いまわすことで極力CUDA API呼び出しを回避 class HashSIFTImpl : public HashSIFT { public: ... private: GpuMat image_, keypoints_, descriptors_, d_bMatrix_; DeviceBuffer bufResponses_; MatmulAndSign matmulAndSign_; }; Copyright © Fixstars Group GPUで使用されるリソース 40

42.

高速化結果 ⚫ 計測環境 CPU Core i9-10900(2.80GHz) 10Core/20T GPU GeForce RTX 3060 Ti および Jetson Xavier 入力画像 SceauxCastle[12]に含まれる11枚の画像 特徴点数 10,000点に固定 計測方法 各画像について100回実行の平均処理時間を計測し、最後に11枚の平均をとる ⚫ 計測結果 実装 HashSIFT256 [msec] HashSIFT512 [msec] リファレンスCPU Single Thread 604.9 740.1 リファレンスCPU Multi Thread 174.6 315.6 0.89 0.99 5.7 6.4 GPU GeForce RTX 3060 Ti GPU Jetson Xavier Copyright © Fixstars Group 41

43.

結果と考察 ⚫ CPUからの高速化倍率が高い ○ CPU Multi Thread → RTX 3060 Ti で200~300倍 ○ SIFT計算はパッチやヒストグラムへの処理など、メモリアクセスがボトルネックと考えられ、 今回それらの処理が全てShared Memory上で完結したのが大きいのではないか ⚫ CUDA化後はHashSIFT256とHashSIFT512で、あまり処理時間の差が無い ○ 処理時間の内訳(下図)を見ると、SIFT計算の割合が大きい ○ SIFT計算はHashSIFT256,512で共通なため、特徴量サイズの違いによる影響が小さい SIFT計算 行列積+2値化 NVIDIA Nsight Systemsで取得したカーネル関数実行のタイムライン Copyright © Fixstars Group 42

44.

おわりに ⚫ ソースコードを公開しています ○ github.com/fixstars/cuda-efficient-features ○ 特徴点検出~特徴量記述までの一連の処理をCUDA実装 ○ ブログ記事:https://proc-cpuinfo.fixstars.com/2023/12/cuda-efficient-features/ ⚫ 今後の展望 ○ SLAMやSfMなど、 cuda-efficient-featuresを応用したアプリケーションの開発、 または既存ライブラリへの移植 Copyright © Fixstars Group 43

45.

参考文献 [1] Rublee, Ethan, et al. "ORB: An efficient alternative to SIFT or SURF." 2011 International conference on computer vision. Ieee, 2011. [2] Lowe, David G. "Distinctive image features from scale-invariant keypoints." International journal of computer vision 60 (2004): 91-110. [3] Suarez, Iago, Jose M. Buenaposada, and Luis Baumela. "Revisiting binary local image description for resource limited devices." IEEE Robotics and Automation Letters 6.4 (2021): 8317-8324. [4] Suárez, Iago, et al. "BEBLID: Boosted efficient binary local image descriptor." Pattern recognition letters 133 (2020): 366-372. [5] https://github.com/iago-suarez/efficient-descriptors [6] Balntas, Vassileios, et al. "Learning local feature descriptors with triplets and shallow convolutional neural networks." Bmvc. Vol. 1. No. 2. 2016. (Triplet Ranking Loss, anchor swap) [7] Mishchuk, Anastasiia, et al. "Working hard to know your neighbor's margins: Local descriptor learning loss." Advances in neural information processing systems 30 (2017). (Hard Negative Mining) [8] Schonberger, Johannes L., et al. "Comparative evaluation of hand-crafted and learned local features." Proceedings of the IEEE conference on computer vision and pattern recognition. 2017. (ETH benchmark) [9] Brown, Matthew, Gang Hua, and Simon Winder. “Discriminative learning of local image descriptors.” IEEE transactions on pattern analysis and machine intelligence 33.1 (2010): 43-57. (Liberty dataset) [10] Mikolajczyk, Krystian, and Cordelia Schmid. "A performance evaluation of local descriptors." IEEE transactions on pattern analysis and machine intelligence 27.10 (2005): 1615-1630. (Oxford dataset) [11] Balntas, Vassileios, et al. "HPatches: A benchmark and evaluation of handcrafted and learned local descriptors." Proceedings of the IEEE conference on computer vision and pattern recognition. 2017. [12] https://github.com/openMVG/ImageDataset_SceauxCastle [13] https://developer.nvidia.com/blog/using-cuda-warp-level-primitives/ Copyright © Fixstars Group 44

46.

補足スライド Copyright © Fixstars Group 45

47.

BADとHashSIFTの性能:Hpatchesベンチマーク[9] ⚫ 以下の3種類の指標により、画像パッチの記述性能を評価 ○ patch verification:パッチの正例と負例の分離性能 ○ image matching:2枚の画像間のマッチング性能 ○ patch retrieval :類似パッチの検索性能 ⚫ ベース手法に対し優位な性能 Copyright © Fixstars Group 46

48.

BADとHashSIFTの性能:ETHベンチマーク[6] ⚫ Madrid Metropolis(1,344画像)をSfMにより3次元復元 ⚫ 復元性能を特徴量毎に評価 ○ 論文ではRegistered(登録画像枚数)とSparse Pts(復元点群数)を重要としている ⚫ 比較結果 ○ BADはORBに対しSparse Ptsが40%程度向上 ○ HashSIFTは軽量なバイナリ特徴でありながら、CNN系に匹敵する性能 ETHベンチマークによるSfM性能比較 ([3]より引用) ORB系 SIFT系 CNN系 Copyright © Fixstars Group 47

49.

Thank you! お問い合わせ窓口 : [email protected] Copyright © Fixstars Group 48