---
title: コンピュータビジョンセミナー vol.4～Visual SLAM / SfMで行われる局所特徴量計算のCUDA高速化～（2024/3/27）
tags:  #cuda #機械学習 #deeplearning #深層学習 #visualslam #sfm #画像解析 #画像処理 #cuda高速化 #自動運転 #コンピュータビジョン #コンピュータビジョンシリーズ  
author: [株式会社フィックスターズ](https://www.docswell.com/user/fixstars)
site: [Docswell](https://www.docswell.com/)
thumbnail: https://bcdn.docswell.com/page/L73RPGX6J5.jpg?width=480
description: 自動移動ロボットに欠かせない 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
published: March 28, 24
canonical: https://www.docswell.com/s/fixstars/ZM1DJ8-20240327
---
# Page. 1

![Page Image](https://bcdn.docswell.com/page/L73RPGX6J5.jpg)

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


# Page. 2

![Page Image](https://bcdn.docswell.com/page/87DYL12M7G.jpg)

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


# Page. 3

![Page Image](https://bcdn.docswell.com/page/VJP62Q1QE8.jpg)

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


# Page. 4

![Page Image](https://bcdn.docswell.com/page/2EV3D5YP7Q.jpg)

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


# Page. 5

![Page Image](https://bcdn.docswell.com/page/57G3Y69QJL.jpg)

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


# Page. 6

![Page Image](https://bcdn.docswell.com/page/4EQ5M4WWJP.jpg)

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


# Page. 7

![Page Image](https://bcdn.docswell.com/page/KJ49VQK171.jpg)

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


# Page. 8

![Page Image](https://bcdn.docswell.com/page/LE1K3PX5JG.jpg)

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


# Page. 9

![Page Image](https://bcdn.docswell.com/page/GEWDP2NWJ2.jpg)

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


# Page. 10

![Page Image](https://bcdn.docswell.com/page/47ZGKD32E3.jpg)

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


# Page. 11

![Page Image](https://bcdn.docswell.com/page/YJ6YQ1Z67V.jpg)

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


# Page. 12

![Page Image](https://bcdn.docswell.com/page/GJ5GR4W5J4.jpg)

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


# Page. 13

![Page Image](https://bcdn.docswell.com/page/9E2YDVQGER.jpg)

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


# Page. 14

![Page Image](https://bcdn.docswell.com/page/D7YG23WG7M.jpg)

代表的な局所特徴量計算手法
⚫ 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


# Page. 15

![Page Image](https://bcdn.docswell.com/page/VENG8V98J8.jpg)

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


# Page. 16

![Page Image](https://bcdn.docswell.com/page/Y79WK62X73.jpg)

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


# Page. 17

![Page Image](https://bcdn.docswell.com/page/G78NRZ59JD.jpg)

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


# Page. 18

![Page Image](https://bcdn.docswell.com/page/L7LPKDYMJR.jpg)

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


# Page. 19

![Page Image](https://bcdn.docswell.com/page/4EM42PN6JW.jpg)

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


# Page. 20

![Page Image](https://bcdn.docswell.com/page/PERWM1DNJ9.jpg)

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


# Page. 21

![Page Image](https://bcdn.docswell.com/page/P7XPVY18JX.jpg)

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


# Page. 22

![Page Image](https://bcdn.docswell.com/page/37KDQ32LED.jpg)

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


# Page. 23

![Page Image](https://bcdn.docswell.com/page/LJ3RPG46E5.jpg)

線形変換+2値化
⚫ SIFT特徴量を線形変換でD次元のベクトルに変換
○
変換行列の係数は事前に学習(後述)
○
n個の特徴点全てのSIFT特徴量をまず計算し、1つの行列に結合後、行列-行列積を行う
⚫ 各要素を閾値処理によって、1/0のバイナリ特徴量に変換
φ
ステップ関数
φ x =ቊ
0 x≤0
1 x&gt;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


# Page. 24

![Page Image](https://bcdn.docswell.com/page/8JDYL1QMJG.jpg)

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


# Page. 25

![Page Image](https://bcdn.docswell.com/page/VEP62QLQ78.jpg)

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


# Page. 26

![Page Image](https://bcdn.docswell.com/page/27V3D5QPJQ.jpg)

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


# Page. 27

![Page Image](https://bcdn.docswell.com/page/5JG3Y6WQEL.jpg)

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


# Page. 28

![Page Image](https://bcdn.docswell.com/page/47Q5M43WEP.jpg)

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


# Page. 29

![Page Image](https://bcdn.docswell.com/page/KE49VQ11J1.jpg)

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


# Page. 30

![Page Image](https://bcdn.docswell.com/page/L71K3PG5EG.jpg)

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


# Page. 31

![Page Image](https://bcdn.docswell.com/page/G7WDP2KWE2.jpg)

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


# Page. 32

![Page Image](https://bcdn.docswell.com/page/4JZGKDZ273.jpg)

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


# Page. 33

![Page Image](https://bcdn.docswell.com/page/YE6YQ186JV.jpg)

実装例：カーネル関数およびカーネル関数呼び出し
__global__ void computeSIFTKernel(const PtrStepSzb image,
const PtrStepSz&lt;KeyPoint&gt; 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&amp; image, const GpuMat&amp; keypoints,
GpuMat&amp; responses)
{
const dim3 block(32, 4); // 32×4のスレッドを生成
const dim3 grid(keypoints.rows, 1); // 特徴点数のブロックを生成
computeSIFTKernel&lt;&lt;&lt;grid, block&gt;&gt;&gt;(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 &lt; DESCRIPTOR_SIZE; k += blockSize1D())
responses(i, k) = descriptor[k];
}
SIFT計算のカーネル関数
※実際の実装より簡略化しています
Copyright © Fixstars Group
32


# Page. 34

![Page Image](https://bcdn.docswell.com/page/GE5GR495E4.jpg)

実装例：局所画像パッチ作成～輝度勾配計算
__device__ void rectifyPatch(const PtrStepSzb image,
const Keypoint&amp; keypoint, uchar patch[PATCH_H][PATCH_W])
{
// キーポイントのスケールとオリエンテーションからアフィン変換行列を作成
const Mat23f M = getAffineTransform(keypoint);
for (int y = threadIdx.y; y &lt; PATCH_H; y += blockDim.y) {
for (int x = threadIdx.x; x &lt; 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 &lt; PATCH_H; y += blockDim.y) {
for (int x = threadIdx.x; x &lt; 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インクリメント
このようにループを記述することで、
スレッド数 &lt; 要素数な処理に柔軟に対応できる
輝度勾配計算
今回の場合、32×4スレッドで32×32のパッチを、
右図のような順序でアクセスしていく
※実際の実装より簡略化しています
Copyright © Fixstars Group
33


# Page. 35

![Page Image](https://bcdn.docswell.com/page/972YDV2G7R.jpg)

実装例：ヒストグラムへの投票
// 勾配強度と勾配方向の計算
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(&amp;histogram[ri
atomicAdd(&amp;histogram[ri
atomicAdd(&amp;histogram[ri
atomicAdd(&amp;histogram[ri
atomicAdd(&amp;histogram[ri
atomicAdd(&amp;histogram[ri
atomicAdd(&amp;histogram[ri
atomicAdd(&amp;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


# Page. 36

![Page Image](https://bcdn.docswell.com/page/DJYG23YGJM.jpg)

実装例：特徴ベクトル作成 (1次元化 + 正規化)
__device__ void createFeatureVector(float histogram[…],
float descriptor[DESCRIPTOR_SIZE])
{
// 1次元化
for (int r = threadIdx.y; r &lt; R_BINS; r += blockDim.y)
for (int c = threadIdx.x; c &lt; C_BINS; c += blockDim.x)
for (int k = 0; k &lt; 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 &lt; DESCRIPTOR_SIZE; i += WARP_SIZE)
sumSq += squared(descriptor[i]);
// Warp ShuffleでWarp内の総和をとる
for (int mask = 16; mask &gt; 0; mask /= 2)
sumSq += __shfl_xor_sync(0xffffffff, sumSq, mask);
const float scale = 1.f / sqrtf(sumSq);
for (int i = threadIdx.x; i &lt; DESCRIPTOR_SIZE; i += WARP_SIZE)
descriptor[i] = scale * descriptor[i];
}
Warp Shuffleによる総和の実装例 ([13]より引用)
}
特徴ベクトル作成 (1次元化 + 正規化)
※実際の実装より簡略化しています
Copyright © Fixstars Group
35


# Page. 37

![Page Image](https://bcdn.docswell.com/page/V7NG8VP8E8.jpg)

ヒストグラム作成時のatomicAdd衝突回数削減
⚫ 複数スレッドが同時に同一ビンにアクセス(衝突)すると待ちが発生
○
ナイーブ実装では32×4スレッドのまとまりが上～下に移動
○
スレッド同士で参照位置が近く、周辺ブロックへの投票も含めると衝突頻度が高いはず
⚫ 参照位置が分散するようにアクセスパターンを工夫
○
効果は微改善程度
for (int y = threadIdx.y; y &lt; PATCH_H; y += blockDim.y) {
for (int x = threadIdx.x; x &lt; 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 &lt; 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


# Page. 38

![Page Image](https://bcdn.docswell.com/page/YJ9WK63XJ3.jpg)

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


# Page. 39

![Page Image](https://bcdn.docswell.com/page/GJ8NRZM9ED.jpg)

行列積
⚫ cuBLASライブラリを使用
void linearProjection(const GpuMat&amp; F, const GpuMat&amp; B, GpuMat&amp; 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(&amp;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,
&amp;alpha,
B.ptr&lt;float&gt;(), B.cols,
F.ptr&lt;float&gt;(), F.cols,
&amp;beta,
H.ptr&lt;float&gt;(), H.step.cols);
// cuBLASハンドルの削除
cublasDestroy_v2(handle);
}
Copyright © Fixstars Group
38


# Page. 40

![Page Image](https://bcdn.docswell.com/page/LJLPKD3MER.jpg)

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++
&gt;
&gt;
&gt;
&gt;
&gt;
&gt;
&gt;
&gt;
0)
0)
0)
0)
0)
0)
0)
0)
&lt;&lt;
&lt;&lt;
&lt;&lt;
&lt;&lt;
&lt;&lt;
&lt;&lt;
&lt;&lt;
&lt;&lt;
7;
6;
5;
4;
3;
2;
1;
0;
dst(i, p) = byte;
}
Copyright © Fixstars Group
39


# Page. 41

![Page Image](https://bcdn.docswell.com/page/47M42PD6EW.jpg)

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


# Page. 42

![Page Image](https://bcdn.docswell.com/page/P7RWM14NE9.jpg)

高速化結果
⚫ 計測環境
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


# Page. 43

![Page Image](https://bcdn.docswell.com/page/PJXPVY28EX.jpg)

結果と考察
⚫ 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


# Page. 44

![Page Image](https://bcdn.docswell.com/page/3JKDQ3ML7D.jpg)

おわりに
⚫ ソースコードを公開しています
○
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


# Page. 45

![Page Image](https://bcdn.docswell.com/page/LE3RPGY675.jpg)

参考文献
[1] Rublee, Ethan, et al. &quot;ORB: An efficient alternative to SIFT or SURF.&quot; 2011 International conference on computer vision. Ieee, 2011.
[2] Lowe, David G. &quot;Distinctive image features from scale-invariant keypoints.&quot; International journal of computer vision 60 (2004): 91-110.
[3] Suarez, Iago, Jose M. Buenaposada, and Luis Baumela. &quot;Revisiting binary local image description for resource limited devices.&quot; IEEE Robotics and
Automation Letters 6.4 (2021): 8317-8324.
[4] Suárez, Iago, et al. &quot;BEBLID: Boosted efficient binary local image descriptor.&quot; Pattern recognition letters 133 (2020): 366-372.
[5] https://github.com/iago-suarez/efficient-descriptors
[6] Balntas, Vassileios, et al. &quot;Learning local feature descriptors with triplets and shallow convolutional neural networks.&quot; Bmvc. Vol. 1. No. 2. 2016. (Triplet
Ranking Loss, anchor swap)
[7] Mishchuk, Anastasiia, et al. &quot;Working hard to know your neighbor&#039;s margins: Local descriptor learning loss.&quot; Advances in neural information processing
systems 30 (2017). (Hard Negative Mining)
[8] Schonberger, Johannes L., et al. &quot;Comparative evaluation of hand-crafted and learned local features.&quot; 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. &quot;A performance evaluation of local descriptors.&quot; IEEE transactions on pattern analysis and machine
intelligence 27.10 (2005): 1615-1630. (Oxford dataset)
[11] Balntas, Vassileios, et al. &quot;HPatches: A benchmark and evaluation of handcrafted and learned local descriptors.&quot; 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


# Page. 46

![Page Image](https://bcdn.docswell.com/page/8EDYL15MEG.jpg)

補足スライド
Copyright © Fixstars Group
45


# Page. 47

![Page Image](https://bcdn.docswell.com/page/V7P62QGQJ8.jpg)

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


# Page. 48

![Page Image](https://bcdn.docswell.com/page/2JV3D56PEQ.jpg)

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


# Page. 49

![Page Image](https://bcdn.docswell.com/page/5EG3Y6NQ7L.jpg)

Thank you!
お問い合わせ窓口 : hr-seminar@fixstars.com
Copyright © Fixstars Group
48


