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

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

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

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

• https://www.docswell.com/s/fixstars/ZRXQ72-20220805
• https://www.docswell.com/s/fixstars/5ENE7J-20221220
• https://www.docswell.com/s/fixstars/ZQ8447-20231027

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

線形変換+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

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

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

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

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

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

• https://www.docswell.com/s/fixstars/K24MYM-20220527

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

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

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

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

実装例：カーネル関数およびカーネル関数呼び出し
__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

実装例：局所画像パッチ作成～輝度勾配計算
__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

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

実装例：特徴ベクトル作成 (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

ヒストグラム作成時の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

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

行列積
⚫ 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

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

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

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

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

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

• https://proc-cpuinfo.fixstars.com/2023/12/cuda-efficient-features/

参考文献 [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

• https://github.com/iago-suarez/efficient-descriptors
• https://github.com/openMVG/ImageDataset_SceauxCastle

補足スライド Copyright © Fixstars Group 45

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

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

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