～リアルタイム画像処理への挑戦～
名古屋大学大学院情報科学研究科
出口 大輔
発表の流れ
GPGPUを始める前に
GPGPUの基礎知識
CUDAの使い方
CUDAを使う前に
プログラミングの予備知識
CUDAを使って“Hello World”
GPGPUにチャレンジ
行列積の計算
テンプレートマッチング
ガウシアンフィルタ
SIFT特徴量の計算
GPGPUって？
GPGPUは何の略？
G
eneral-
P
urpose computation on
GPU
s
GPUを汎用計算に利用しようという試み
現在は「GPUコンピューティング」とも呼ばれる
0 200 400 600 800 1000 1200 1400 1600 1900 1900 1900 1900 1900 1900 1900 1900 1900
CPUとGPUの性能比較
2003
2004
2005
2006
2008
※NVIDIA CUDA Programming Guide 4.0より引用
[GFLOP/s]
ピ
ー
ク
性
能
NVIDIA GPU
Intel
CPU
Geforece GTX280
Geforece 8800GTX
Quad Core
Xeon 3.2GHz
Core2 Duo
3.0GHz
2007
2009
2010
Geforece GTX480
GPGPUって？
GPGPUは何の略？
G
eneral-
P
urpose computation on
GPU
s
GPUを汎用計算に利用しようという試み
現在は「GPUコンピューティング」とも呼ばれる
どうしてGPGPUが注目されているのか？
GPUの計算性能が飛躍的に向上
最新のGPUは 1.5 TFLOPS 以上の演算性能（CPUの約10倍）
GPUはCPUと比較して並列計算に優れている
GeForce GTX580 では 512 コアによる並列計算が可能
手頃な価格で入手可能
GPGPUの活用例
動画像処理
CyberLink PowerDirector 7
ビデオエフェクトのレンダリングを高速化
TMPGEnc 4.0 XPress
フィルター処理・デコード処理の高速化
画像処理
Adobe Photoshop CS4 / CS5，Adobe Premire CS5
各種フィルタ処理の高速化
OpenCV
各種画像処理の高速化
数値計算
MATLAB
FFTの高速化
GPUの歴史: 1999年以前
1970年 ～ 1990年
初期開発の時代
ソフトウェアによるグラフィックス処理
プログラム可能なGPUに関する初期研究
Ikonas System [1], Pixel Machine [2]
1990年 ～ 1999年
GPU技術の黎明期
3Dグラフィックス・アクセラレータの開発
グラフィックス向けライブラリの開発
OpenGL（1992年～）, DirectX（1995年～）
[1] J. N. England, “A system for interactive modeling of physical curved surface
objects,” Proceedings of SIGGRAPH 78, pp.336-340. 1978
[2] M. Potmesil and E. M. Hoffert, “The Pixel Machine: A Parallel Image
Computer,” Proceedings of SIGGRAPH89, pp.69-78, 1989
GPUの歴史: 1999年～
“GPU”の誕生
NVIDIA社の GeForce 256 の登場
ハードウェア T＆L
※
をサポート
CPU負荷を大幅に削減
グラフィックスパイプライン
ハードウェア固定の処理
自由表現な表現は不可
頂点処理
（ハードウェア固定）
ジオメトリ処理
（クリッピング等）
ラスタライズ処理
画面出力
※Transform & Lighting
ピクセル処理
GPUの歴史: 2003年～
プログラマブルシェーダの登場
Vertex Shader
頂点座標の変換処理
Pixel Shader
画素の輝度計算処理
シェーダ言語の進化
アセンブラから高級言語へ
Cg, HLSL, GLSL
柔軟な映像表現が可能に
グラフィックス以外への応用
GPGPUへの関心が高まる
Vertex Shader
（プログラマブル）
ジオメトリ処理
（クリッピング等）
ラスタライズ処理
画面出力
Pixel Shader
（プログラマブル）
GPUの歴史: 2007年～
GPGPUの開発環境の整備
CUDA (NVIDIA)
ATI Stream (AMD)
OpenCL
DirectCompute
GPGPU時代の到来
物理シミュレーション
数値計算
信号解析
画像処理・認識
･･･
Vertex Shader
（プログラマブル）
Geometry Shader
（プログラマブル）
ラスタライズ処理
画面出力
Pixel Shader
（プログラマブル）
統合型シェーダ
GPUの歴史: 2009年～
Fermiアーキテクチャの登場
汎用計算向けのアーキテクチャ
※
L1/L2キャッシュの搭載
複数カーネルの同時実行のサポート
倍精度浮動小数点演算の高速化
ECCメモリのサポート
アトミックなメモリ操作の高速化
C++のフルサポート
GPGPU関連のライブラリ&ツールの充実
CUBLAS, CUFFT, Thrust, NPP
Parallel Nsight, Visual Profiler
CUDAって何？
CUDA（Compute Unified Device Architecture）
発音は「クーダ」もしくは「キューダ」
NVIDIA社が提供するGPUを利用する
ための統合開発環境
GeForce
8以降のハードウェアで利用可能
グラフィックス処理APIの知識は不要
CPUでプログラムを実行する感覚
C/C++を用いてプログラムの開発が可能
既存のアルゴリズムの移植も比較的容易
CUDAが利用可能な環境
CUDA対応のグラフィックスカード
OSの対応状況（CUDA 4.0）
Windows XP（32bit, 64bit）
Windows Vista（32bit, 64bit）
Windows Server 2008（32bit, 64bit）
Linux（32bit, 64bit）
Mac OS
GeForce
GTX 580, GTX 480, GTX 260, 8800シリーズ，他
Quadro
Plex 2200 D2, FX 5800, FX 5600, 5000, 6000, 他
Tesla
S2050, C2070, S1070, C1060, S870, D870, C870
※TeslaはHPCに特化した製品であり映像出力を持たない
CUDAを使うための準備（Windows編）
「CUDA ZONE」から次をダウンロード
NVIDIA Driver
CUDA対応のビデオドライバー
CUDA Toolkit 4.0
※
コンパイラ（nvcc）
CUBLASやCUFFT
ライブラリ
ドキュメント
CUDA SDK 4.0
サンプル
※ 2011年5月25日にリリース
http://www.nvidia.com/object/cuda_home.html
CUDAを動かしてみよう！！
Volumetric Particle Shadows
※「CUDA SDK」内のサンプルの実行結果
Image Denoising
CUDA対応のGPU（G80, GT200）
ストリーミング・マルチプロセッサ（SM）
8 個のスカラープロセッサ（SP）
16KBの共有メモリ
スレッド間の同期機構
マルチプロセッサ内でのみ可能
マルチプロセッサ間での同期には
CPU処理が必要
GeForce GTX280 の場合
30 基のマルチプロセッサを搭載
1GByte以上のグローバルメモリ
SP
SP
SP
SP
SP
SP
SP
SP
マルチプロセッサ #1
SP
SP
SP
SP
SP
SP
SP
SP
マルチプロセッサ #2
CUDA対応のGPU（Fermi）
ストリーミング・マルチプロセッサ（SM）
32
個のスカラープロセッサ（SP）
48
KBの共有メモリ
スレッド間の同期機構
マルチプロセッサ内でのみ可能
マルチプロセッサ間での同期には
CPU処理が必要
GeForce GTX480 の場合
15
基のマルチプロセッサを搭載
1GByte以上のグローバルメモリ
SM #1
SP SP
SP
SP
SP
SP
SP SP
SP SP
SP
SP
SP
SP
SP SP
SP SP
SP
SP
SP
SP
SP SP
SP SP
SP
SP
SP
SP
SP SP
SM #2
SP SP
SP
SP
SP
SP
SP SP
SP SP
SP
SP
SP
SP
SP SP
SP SP
SP
SP
SP
SP
SP SP
SP SP
SP
SP
SP
SP
SP SP
CUDAのプログラミングモデル
GPUの特徴
多数のスレッドが高い並列性をもって処理を実行
GPU のみでプログラムを実行できない
CPUとの連携が不可欠
CUDAにおける計算の流れ
GPUを並列演算可能なデバイスとして扱う
複数のスレッドを同時に実行できる外部CPU
階層的にスレッドを管理
スレッドのまとまり = ブロック
ブロックのまとまり = グリッド
問題を分割して計算する際に便利
スレッドを3次元的に配置
各スレッドの ID を X, Y, Z の3要素で表現
スレッドのまとまりをブロックとして管理
グリッド
ブロック(0,1,1)
ブロック(0,0,1)
スレッド
(0,0,1)
スレッド
(0,1,1)
スレッド
(0,2,1)
スレッド
(1,0,1)
スレッド
(1,1,1)
スレッド
(1,2,1)
ブロック(1,0,1)
スレッド
(0,0,1)
スレッド
(0,1,1)
スレッド
(0,2,1)
スレッド
(1,0,1)
スレッド
(1,1,1)
スレッド
(1,2,1)
ブロック(1,1,1)
ブロック(0,1,0)
ブロック(1,1,0)
ブロック(1,0,0)
スレッド
(0,0,0)
スレッド
(0,1,0)
スレッド
(0,2,0)
スレッド
(1,0,0)
スレッド
(1,1,0)
スレッド
(1,2,0)
スレッド
(0,0,0)
スレッド
(1,0,0)
スレッド
(2,0,0)
スレッド
(0,1,0)
スレッド
(1,1,0)
スレッド
(2,1,0)
CUDAにおけるスレッド管理
ブロック(0,0,0)
スレッド
(0,0,0)
スレッド
(0,1,0)
スレッド
(0,2,0)
スレッド
(1,0,0)
スレッド
(1,1,0)
スレッド
(1,2,0)
スレッド
(0,0,0)
スレッド
(1,0,0)
スレッド
(2,0,0)
スレッド
(0,1,0)
スレッド
(1,1,0)
スレッド
(2,1,0)
階層的スレッドの利用方法
一般的な画像処理で利用する場合
ブロック内のスレッド数を決定
16×16 = 256 スレッド
画像内にブロックを配置
ブロック内のスレッド
が各画素を処理
計算範囲の求め方
ブロックID
スレッドID
ブロック内スレッド数
T
00T
10T
N0T
01ブロック#1
T
0MT
NM T00 T10 TN0 T01 T0M TNMスレッド
(0,0,0)
スレッド
(1,0,0)
スレッド
(2,0,0)
スレッド
(0,1,0)
スレッド
(1,1,0)
スレッド
(2,1,0)
ブロック(0,0,0)
CUDAにおける計算の流れ
ブロック(1,0,0)
ブロック(0,1,0)
ブロック(1,1,0)
処理１
グリッド 1
ブロック (0,0,0)
ブロック (1,0,0)
グリッド 2
処理２
スレッド
(0,0,0)
スレッド
(1,0,0)
スレッド
(2,0,0)
スレッド
(0,1,0)
スレッド
(1,1,0)
スレッド
(2,1,0)
CUDAのメモリモデル
共有メモリ
スレッド (0,0,0)
ローカルメモリ
レジスタ
スレッド(1,0,0)
ブロック (0,0,0)
共有メモリ
スレッド(0,0,0)
スレッド(1,0,0)
ブロック(1,0,0)
コンスタントメモリ
テクスチャメモリ
グローバルメモリ
ローカルメモリ
レジスタ
ローカルメモリ
レジスタ
ローカルメモリ
レジスタ
CUDAで利用可能なメモリ
グローバルメモリ
大量のメモリを利用可能
低速なメモリアクセス（400～600クロック必要）
Fermiアーキテクチャではキャッシュ機構（L1/L2）を搭載
共有メモリ
レジスタと同じ速度でアクセス可能
マルチプロセッサ1基あたり16KB（Fermiは最大48KB）
Fermiの場合は一部をL1キャッシュとして利用可能
テクスチャメモリ
キャッシュ機構による高速なアクセスが可能
ハードウェア線形補間や正規化テクスチャ座標が利用可能
コンスタントメモリ
キャッシュ機構による高速なアクセスが可能
マルチプロセッサ1基あたり64KB
各GPUで使用できる機能の違い
使用できる機能を Compute Capability で区別
現在 1.0 ～ 2.1 のGPUが存在
Compute Capability による機能の違い
1.0：初期リリース
GeForce 8800GTX, Quadro FX 5600, 他
1.1：アトミックなメモリ操作のサポート
GeForce 9800 GTX, Quadro FX 3700, 他
1.3：倍精度浮動小数点演算のサポート
GeForce GTX 280, Quadro FX 5800, 他
2.x ：グローバルメモリのキャッシュをサポート
GeForce GTX480，Quadro 5000，他
Compute Capabilityによる機能の違い
Compute Capability
1.0
1.1
1.2
1.3
2.x
3次元グリッド
×
○
倍精度浮動小数
×
○
32ビット整数の
アトミック演算
×
○
64ビット整数の
アトミック演算
×
○
単精度浮動小数の
アトミック演算
×
○
CUDAにおける制限事項
Compute Capability
1.0
1.1
1.2
1.3
2.x
1ブロックあたりの
スレッド数
512
1024
1マルチプロセッサ
あたりのスレッド数
768
1024
1536
1マルチプロセッサ
あたりのレジスタ数
8192
16384
32768
コンスタントメモリ
64KB
共有メモリ
16KB
48KB
2Dテクスチャ
2
16
×2
15
2
16
×2
16
CUDAを使って“Hello World!!”
Hello World!! を表示するサンプル（main.cu）
__global__ void hello( char *data ) {
char *text = "Hello World!!¥n";
data[ threadIdx.x ] = text[ threadIdx.x ]; }
int main( int argc, char *argv[] ) {
char *dData, hData[ 15 ];
cudaMalloc( ( void ** )&dData, sizeof( char ) * 15 ); dim3 nThreads( 15, 1 );
dim3 nBlocks( 1, 1 );
hello<<< nBlocks, nThreads >>>( dData );
cudaMemcpy( hData, dData, sizeof( char ) * 15, cudaMemcpyDeviceToHost ); printf( "%s", hData );
cudaFree( dData ); return( 0 );
CUDAを使って“Hello World!!”
Hello World!! を表示するサンプル（main.cu）
__global__ void hello( char *data ) {
char *text = "Hello World!!¥n";
data[ threadIdx.x ] = text[ threadIdx.x ]; }
int main( int argc, char *argv[] ) {
char *dData, hData[ 15 ];
cudaMalloc( ( void ** )&dData, sizeof( char ) * 15 ); dim3 nThreads( 15, 1 );
dim3 nBlocks( 1, 1 );
hello<<< nBlocks, nThreads >>>( dData );
cudaMemcpy( hData, dData, sizeof( char ) * 15, cudaMemcpyDeviceToHost ); printf( "%s", hData ); cudaFree( dData ); return( 0 ); }
_
_global_
_, threadIdx
dim3, <<<…>>>
CUDAで拡張された部分
CUDAを使って“Hello World!!”
Hello World!! を表示するサンプル（main.cu）
__global__ void hello( char *data ) {
char *text = "Hello World!!¥n";
data[ threadIdx.x ] = text[ threadIdx.x ]; }
int main( int argc, char *argv[] ) {
char *dData, hData[ 15 ];
cudaMalloc( ( void ** )&dData, sizeof( char ) * 15 ); dim3 nThreads( 15, 1 );
dim3 nBlocks( 1, 1 );
hello<<< nBlocks, nThreads >>>( dData );
cudaMemcpy( hData, dData, sizeof( char ) * 15, cudaMemcpyDeviceToHost ); printf( "%s", hData );
cudaFree( dData ); return( 0 );
CUDAの言語拡張（1）
言語拡張により追加された修飾子
関数に対する修飾子
変数に対する修飾子
_ _global_ _
CPU から呼び出され，GPU で実行される関数
_ _device_ _
GPU から呼び出され，GPU で実行される関数
_ _host_ _
CPU から呼び出され，CPU で実行される関数
_ _constant_ _ GPU 上のコンスタントメモリに存在する変数
_ _shared_ _
GPU 上の共有メモリに存在する変数
_ _device_ _
GPU 上のメモリに存在する変数
CUDAを使って“Hello World!!”
Hello World!! を表示するサンプル（main.cu）
__global__ void hello( char *data ) {
char *text = "Hello World!!¥n";
data[ threadIdx.x ] = text[ threadIdx.x ]; }
int main( int argc, char *argv[] ) {
char *dData, hData[ 15 ];
cudaMalloc( ( void ** )&dData, sizeof( char ) * 15 ); dim3 nThreads( 15, 1 );
dim3 nBlocks( 1, 1 );
hello<<< nBlocks, nThreads >>>( dData );
cudaMemcpy( hData, dData, sizeof( char ) * 15, cudaMemcpyDeviceToHost ); printf( "%s", hData );
cudaFree( dData ); return( 0 );
CUDAの言語拡張（2）
CUDAで利用可能な組み込み型
GPU内のスレッドを識別する組み込み変数
dim3
整数 x, y, z からなる 3 次元ベクトル
（スレッド数やブロック数の指定に利用）
uchar2, int2, float2, …
x, y からなる 2 次元ベクトル
uchar3, int3, float3, …
x, y, z からなる 3 次元ベクトル
uchar4, int4, float4, …
x, y, z, w からなる 4 次元ベクトル
gridDim
グリッドの次数
blockIdx
スレッドが属するブロックのインデックス
blockDim
スレッドが属するブロックの次数
gridDim
= 2×2×2
スレッドを識別する組み込み変数
グリッド
ブロック(0,1,1)
ブロック(1,0,1)
スレッド
(0,0,1)
スレッド
(0,1,1)
スレッド
(0,2,1)
スレッド
(1,0,1)
スレッド
(1,1,1)
スレッド
(1,2,1)
ブロック(1,1,1)
ブロック(0,1,0)
ブロック(1,1,0)
ブロック(1,0,0)
スレッド
(0,0,0)
スレッド
(0,1,0)
スレッド
(0,2,0)
スレッド
(1,0,0)
スレッド
(1,1,0)
スレッド
(1,2,0)
スレッド
(0,0,0)
スレッド
(1,0,0)
スレッド
(2,0,0)
スレッド
(0,1,0)
スレッド
(1,1,0)
スレッド
(2,1,0)
ブロック(0,0,1)
スレッド
(0,0,1)
スレッド
(0,1,1)
スレッド
(0,2,1)
スレッド
(1,0,1)
スレッド
(1,1,1)
スレッド
(1,2,1)
ブロック
(0,0,0)
blockIdx
blockDim
= 3×2×2
スレッド
(0,0,0)
スレッド
(0,1,0)
スレッド
(0,2,0)
スレッド
(1,0,0)
スレッド
(1,1,0)
スレッド
(1,2,0)
スレッド
(1,0,0)
スレッド
(2,0,0)
スレッド
(0,1,0)
スレッド
(1,1,0)
スレッド
(2,1,0)
スレッド
(0,0,0)
threadIdx
C/C++プログラミングとの相違点
CPUとGPUで実行するコードを明確に区別
CPUで実行する場合は “_ _host_ _”（省略可）を付与
GPUで実行する場合は “_ _global_ _” を付与
CPUとGPUで処理可能なメモリ空間の違い
CPUとGPU間でメモリ転送が必要
スレッド数を指定してGPU上の関数を呼び出し
関数名<<< ブロック数, スレッド数 >>>( … );
ブロック数とスレッド数は実行時に指定可能
問題サイズに合わせて調整可能
CUDAを使って“Hello World!!”
Hello World!! を表示するサンプル（main.cu）
__global__ void hello( char *data ) {
char *text = "Hello World!!¥n";
data[ threadIdx.x ] = text[ threadIdx.x ]; }
int main( int argc, char *argv[] ) {
char *dData, hData[ 15 ];
cudaMalloc( ( void ** )&dData, sizeof( char ) * 15 ); dim3 nThreads( 15, 1 );
dim3 nBlocks( 1, 1 );
hello<<< nBlocks, nThreads >>>( dData );
cudaMemcpy( hData, dData, sizeof( char ) * 15, cudaMemcpyDeviceToHost ); printf( "%s", hData ); cudaFree( dData ); return( 0 ); }
GPU上で実行される関数
GPU上のメモリを解放
CPUからGPUへメモリ転送
GPU上にメモリを確保
ブロック内に配置する
スレッド数 15×1 = 15
総ブロック数 1×1 = 1
スレッドのレジスタ数を調査
Visual Studio 2008 コマンドプロンプトを起動
コマンドプロンプト上で main.cu をコンパイル
合計レジスタ数をチェック
C:¥Your¥Path>nvcc main.cu --ptxas-options=-v --compile
main.cu
ptxas info : Compiling entry function '_Z5helloPc' for 'sm_10'
ptxas info :
Used 2 registers
, 8+16 bytes smem, 15 bytes cmem[0]
レジスタ数
2 × 15 = 30
Thrustライブラリって？
CUDAで利用できるテンプレートライブラリ
CUDAとOpenMPに対応したC++ STLの並列版
URL: http://thrust.googlecode.com/
Thrustライブラリの特徴
コンテナ
CPUとGPUのデータをコンテナとして管理
アルゴリズム
コンテナに対してアルゴリズムを適用可能
直感的なインターフェース
CUDAの複雑なAPIに関する知識は不要
メモリの確保／解放／転送がとても簡単
Thrustライブラリの使い方
GPU上で「1,2,3,…,100」の総和を計算
#include <thrust/host_vector.h> #include <thrust/device_vector.h> #include <thrust/sequence.h> #include <thrust/reduce.h>int main( int argc, char *argv[] ) {
// CPU上のメモリを確保
thrust::host_vector< int > hvec( 100 );
// データを初期化 1, 2, 3, ...
thrust::sequence( hvec.begin( ), hvec.end( ), 1 );
// GPU上のメモリを確保
thrust::device_vector< int > dvec( 100 );
// CPU -> GPU のメモリ転送
dvec = hvec;
int val = thrust::reduce( dvec.begin( ), dvec.end( ), 0, thrust::plus< int >( ) );
printf( "%d¥n", val );
return( 0 ); }
コンテナの使い方（1）
メモリの確保
thrust::host_vector< int > hvec( 20 );
CPU上のメモリに要素数20の「int」配列を確保
thrust::device_vector< float > dvec( 100 );
GPU上のメモリに要素数100の「float」配列を確保
メモリの解放
コンテナオブジェクトの消滅時に自動解放
CPUの場合：free( )
GPUの場合：cudaFree( )
明示的なメモリ解放：dvec.clear( )
コンテナの使い方（2）
メモリの転送
代入操作でメモリ転送が可能
コンテナの各要素へのアクセス
配列と同様にアクセス可能
thrust::host_vector< float >
hvec( 20 );
thrust::device_vector< float > dvec( 20 );
dvec = hvec;
hvec[ 5 ] = 100.0f;
dvec[ 1 ] = 23.4f;
イテレータ（1）
コンテナの要素を指すポインタ
STLのイテレータと同じように使用可能
thrust::device_vector< int > dvec( 4 );
thrust::device_vector< int >::iterator ite = dvec.begin( );
0
1
2
3
dvec.begin( )
dvec.end( )
dvec.begin( ) + 3
イテレータ（2）
コンテナの各要素へのアクセス
thrust::device_vector< int > dvec( 4 );
thrust::device_vector< int >::iterator ite = dvec.begin( );
*ite = 10;
++ite;
*ite = 25;
dvec[ 0 ] = 10 と等価
アルゴリズムと使用方法
CPUとGPUのデータに対して同じ様に適用可能
GPUの場合は
GPUを使用して並列処理
される
アルゴリズムの使用例
コンテナのデータすべてに1を代入する
thrust::fill( dvec.begin( ), dvec.end( ), 1 );
コンテナのデータに 1, 2, 3, … の値で初期化する
thrust::sequence( dvec.begin( ), dvec.end( ) );
1 1 1 1
1 1 1 1
Parallel Reduction（1）
コンテナ内の要素を1つの値に集約する処理
GPUのスレッドが並列に処理を実行
例：総和計算の場合
36
1 2 3 4 5 6 7 8
GPU上のデータ
3
7
11
15
10
26
36
CUDAスレッドの同期
CUDAスレッドの同期
Parallel Reduction（2）
Thrustライブラリを用いた集約処理
総和の計算
int val = thrust::reduce( dvec.begin( ), dvec.end( ) );
最大値の計算
int val = thrust::reduce( dvec.begin( ), dvec.end( ),
0, thrust::maximum< int >( ) );
1 2 3 4 5 6 7
10
55
CUDAとThrustを組み合わせる
Hello World!! を表示するサンプル
#include <thrust/host_vector.h> #include <thrust/device_vector.h> __global__ void hello( char *data ) {
char *text = "Hello World!!¥n";
data[ threadIdx.x ] = text[ threadIdx.x ]; }
int main( int argc, char *argv[] ) {
thrust::host_vector< char > hvec( 15 ); thrust::device_vector< char > dvec( 15 ); dim3 nThreads( 15, 1 );
dim3 nBlocks( 1, 1 );
hello<<< nBlocks, nThreads >>>( thrust::raw_pointer_cast( &dvec[ 0 ] ) ); hvec = dvec; printf( "%s", &hvec[ 0 ] ); return( 0 ); }
CPUからGPUへメモリ転送
GPU上のメモリを指す
ポインタを取得
GPGPUにチャレンジ
行列積の計算
CUDAにおける基本的な実装方法
共有メモリ
テンプレートマッチング
2次元テクスチャメモリ
ハードウェア線形補間
正規化テクスチャ座標
Parallel Reduction アルゴリズム
ガウシアンフィルタ
1次元テクスチャメモリ
コンスタントメモリ
SIFT特徴量の計算
1・2次元テクスチャメモリ
コンスタントメモリ
GPUによる行列積の計算
行列積計算の特徴
各要素はすべて独立に計算可能
GPUによる並列計算が可能
同じメモリ領域へ頻繁にアクセス
メモリアクセス速度がボトルネック
GPU上での実装方法
1.
行列積の計算式に忠実な実装（GPU#1）
2.
共有メモリを利用した高速化（GPU#2）
ブロック単位（部分行列）で行列積を計算
共有メモリをキャッシュとして利用
行列積C＝A×Bの流れ
1.
行列A, B, CのメモリをGPU上に確保
2.
行列A, BのデータをGPUに転送
3.
行列Cの各要素c
mn
を計算
N 次正方行列を仮定
ただし，N は16の倍数
4.
計算結果をCPUに転送
N
k
kn
mk
mn
a
b
c
1
行列A
行列C
行列B
mn
c
mk
a
kn
b
m行
n列
N
N
処理の流れ
CPU側の処理
行列の初期化
GPU側の処理
行列Cの各要素を計算
行
列
の
初
期
化
CPU
メ
モ
リ
確
保
GPU
行
列
積
の
計
算
GPU
デ
ー
タ
転
送
CPU GPU
デ
ー
タ
転
送
GPU CPU
int main( int argc, char *argv[] )
{
int N = 512; // 行列A, B, Cのサイズ
float *hA, *hB, *hC; // CPU（host）側で利用するメモリへのポインタ
float *dA, *dB, *dC; // GPU（device）側で利用するメモリへのポインタ
/* CPU側のメモリを確保*/
/* GPU側のメモリを確保*/
/* CPU側のメモリをGPU側へ転送 */
/* 実行するGPUのスレッド数，ブロック数を設定 */
/* GPUのカーネルを実行し，C=A×B の結果を dC に格納 */
/* GPUの計算結果をCPU側へ転送 */
/* CPUとGPUそれぞれのメモリを解放 */
return( 0 );
}
行列積計算のCPU処理（ソース）
行列積計算のCPU処理（1）
行列積の計算で使用する変数
CPUとGPUでメモリを確保
CPU：malloc, new, cudaMallocHost でメモリを確保
GPU：cudaMalloc によりグローバルメモリを確保
int N = 512; // 行列A, B, Cのサイズ
float *hA, *hB, *hC; // CPU（host）側で利用するメモリへのポインタ float *dA, *dB, *dC; // GPU（device）側で利用するメモリへのポインタ
/* CPU側のメモリを確保 */
hA = ( float * )malloc( N * N * sizeof( float ) ); hB = ( float * )malloc( N * N * sizeof( float ) ); hC = ( float * )malloc( N * N * sizeof( float ) ); /* GPU側のメモリを確保 */
cudaMalloc( ( void ** )&dA, N * N * sizeof( float ) ); cudaMalloc( ( void ** )&dB, N * N * sizeof( float ) ); cudaMalloc( ( void ** )&dC, N * N * sizeof( float ) );
行列積計算のCPU処理（2）
CPUとGPU間のメモリ転送
CPU→GPU
cudaMemcpy( dst, src, size, cudaMemcpyHostToDevice );
GPU→CPU
cudaMemcpy( dst, src, size, cudaMemcpyDeviceToHost );
CPUとGPUで確保したメモリの解放
CPU：free, delete, cudaFreeHost でメモリを解放
GPU：cudaFree でグローバルメモリを解放
転送先
転送バイト数
転送元
転送方向
行列積計算のCPU処理（3）
GPU上で処理を実行
計算範囲の設定
N 次正方行列を仮定
ただし N は16の倍数
行列Cをサイズ16×16のブロックに分割
ブロック内の各要素を各スレッドが計算
合計M×M 個のブロックを配置
各スレッドが c
mn
を計算
n = threadIdx.x + blockIdx.x × blockDim.x;
m = threadIdx.y + blockIdx.y × blockDim.y;
GPU上で関数を実行
multiply<<< ブロック数, スレッド数 >>>( … );
行列C
16
16
(threadIdx.y, threadIdx.x)
(blockIdx.y, blockIdx.x)
行方向のブロック数 (=16)
int main( int argc, char *argv[] )
{
int N = 512; // 行列A, B, Cのサイズ
float *hA, *hB, *hC; // CPU（host）側で利用するメモリへのポインタ
float *dA, *dB, *dC; // GPU（device）側で利用するメモリへのポインタ
/* CPU側のメモリを確保*/
/* GPU側のメモリを確保*/
/* CPU側のメモリをGPU側へ転送 */
/* 実行するGPUのスレッド数，ブロック数を設定 */
dim3 nThreads( 16, 16 );
dim3 nBlocks( N / nThreads.x, N / nThreads.y );
/* GPUのカーネルを実行し，C=A×B の結果を dC に格納 */
multiply<<< nBlocks, nThreads >>>( dA, dB, dC, N );
/* GPUの計算結果をCPU側へ転送 */
/* CPUとGPUそれぞれのメモリを解放 */
return( 0 );
}
行列積計算のCPU処理（まとめ）
GPU上で処理を実行
GPU上のメモリを指定
実装方法（GPU#1）
行列積C＝A×Bの各要素c
mn
を計算
要素c
mn
をGPU上で並列計算
__global__ void multiply1( float *dA, float *dB, float *dC, int N )
{
int n = threadIdx.x + blockIdx.x * blockDim.x;
int m = threadIdx.y + blockIdx.y * blockDim.y;
float sum = 0.0f;
for( int k = 0 ; k < N ; k++ )
{
sum += dA[ m + k * N ] * dB[ k + n * N ];
}
dC[ m + n * N ] = sum;
}
N
k
kn
mk
mn
a
b
c
1
行列A
行列C
行列B
mn
c
mk
a
kn
b
m行
N
N
各要素c
mn
を部分行列（ブロック）の積で計算
ブロックサイズが16×16の場合
ブロック内で同じメモリを参照
高速な共有メモリを利用
グローバルメモリへの
アクセスを削減
2×N×16
2
回
◦
2×N×16
回
実装方法（GPU#2）
16
1
)
1
(
16
16
N
t
t
t
k
kn
mk
mn
a
b
c
行列A
行列C
行列B
mn
c
mk
a
kn
b
N
N
16
16
16
16
16分の1
実装方法（GPU#2）
__global__ void multiply2( float *dA, float *dB, float *dC, int N )
{
int n = threadIdx.x + blockIdx.x * blockDim.x;
int m = threadIdx.y + blockIdx.y * blockDim.y;
float sum = 0.0f;
for( int k = 0 ; k < N ; k += 16 )
{
__shared__ float tA[ 16 ][ 16 ];
__shared__ float tB[ 16 ][ 16 ];
tA[ threadIdx.y ][ threadIdx.x ] = dA[ m + ( k + threadIdx.x ) * N ];
tB[ threadIdx.x ][ threadIdx.y ] = dB[ ( k + threadIdx.y ) + n * N ];
__syncthreads( );
for( int t = 0 ; t < 16 ; t++ )
{
sum += tA[ threadIdx.y ][ t ] * tB[ threadIdx.x ][ t ];
}
__syncthreads( );
}
dC[ m + n * N ] = sum;
}
共有メモリの宣言
ブロック内のスレッドを同期
ブロック内のスレッドを同期
各スレッドが独立
して行列積を計算
各スレッドが独立して
メモリアクセス
行列A 行列C 行列B mnc
mka
knb
N
N
16 16 16 16計算時間の比較（1）
CPU, GPU#1, GPU#2 の計算時間を比較
使用計算機
CPU: Intel Core i7-X980（3.33GHz）
OpenMPを使用して6スレッドで並列計算
GPU: NVIDIA Geforce GTX480
マルチプロセッサ数: 15（SP数 = 15×32 = 480）
Memory：1.5 GB
0 50 100 150 200 250 300 350 400 450 500 16 32 48 64 80 96 11 2 12 8 14 4 16 0 17 6 19 2 20 8 22 4 24 0 27 2 28 8 30 4 32 0 33 6 35 2 36 8 40 0 41 6 43 2 44 8 46 4 48 0 49 6 52 8 54 4 56 0 59 2 60 8 62 4 65 6 67 2 68 8 72 0 73 6 75 2 76 8 78 4 80 0
計算時間の比較（2）
行列のサイズ （N
次正方行列）
[msec.]
CPU
GPU#1
GPU#2
N = 560
計算時間（N
= 560 の場合）
･CPU
#1
: 91.7 msec.
･GPU#1: 18.3 msec. (×5.0)
･GPU#2:
0
6.7 msec. (
×13.7
)
まとめ（行列積）
GPU上での実装方法
1.
行列積の計算式に忠実な実装（GPU#1）
2.
共有メモリを利用した高速化（GPU#2）
大きさ16×16の部分行列積に分割
共有メモリをキャッシュとして利用
CPUとGPUそれぞれでの計算速度を比較
CPU ＜ GPU#1 ＜ GPU#2
メモリアクセスの工夫により約34倍の高速化
テンプレートマッチング（1）
画像処理の分野で広く用いられている手法
基板の品質検査
画像中の特定物体（人物など）の検出
入力画像中からテンプレートに類似する部分を探索
入力画像中に窓を設定
窓を移動させながら
類似度を評価
テンプレート
窓
比較
テンプレートマッチング（2）
基本的なテンプレートマッチングの戦略
入力画像中の部分画像とテンプレートの類似度を評価
窓を移動させながら類似度計算
位置（X軸方向, Y軸方向）
拡大／縮小（スケール変化）
回転
類似度の例
SAD（Sum of Absolute Difference）
SSD（Sum of Squared Difference）
NCC（Normalized Cross Correlation）
窓の大きさに比例して計算コストが増加
膨大な回数の類似度評価が必要
テンプレートマッチングの実装方法
テンプレートマッチングの特徴
類似度は各位置で独立に計算可能
GPUによる並列計算が可能
頻繁に画像メモリへアクセス
テクスチャメモリのキャッシュ機能を利用
複数スケールでテンプレートと入力画像を比較
正規化テクスチャ座標とハードウェア線形補間の利用
解像度に依存しないメモリアクセス
･･･
スケール
処理の流れ
CPU側の処理
入力画像とテンプレートの読み込み
GPUスレッドの同期
GPU側の処理
類似度計算
スケールの変更
画
像
読
み
込
み
CPU
メ
モ
リ
確
保
GPU
類
似
度
計
算
GPU
ス
レ
ッ
ド
同
期
CPU
デ
ー
タ
転
送
CPU GPU
デ
ー
タ
転
送
GPU CPU
テクスチャメモリとは？
GPU上の読み取り専用の特殊なメモリ領域
キャッシュを利用した高速なアクセスが可能
2次元アクセスに対して効率的
ハードウェアを利用した高速な線形補間が可能
テクスチャメモリの定義
texture<DataType, Type, ReadMode> texRef;
DataType: データの型（int，float 等）
Type: テクスチャの種類（1次元テクスチャなど）
cudaTextureType1D，cudaTextureType2D，…
ReadMode: 値の範囲
cudaReadModeElementType
: 各データ型の値を使用
cudaReadModeNormalizedFloat
: 0 ～ 1に正規化（符号付きは-1～1）
メモリ確保と転送（1）
入力画像のメモリをGPU上に確保
テクスチャの定義（画素表現：RGBA）
cudaArrayを利用して2次元テクスチャを確保
CPUからGPUへメモリ転送
テクスチャメモリへの対応付け
cudaArray *imgArray;
cudaChannelFormatDesc c1 = cudaCreateChannelDesc< uchar4 >( );
cudaMallocArray( &imgArray, &c1, width, height );
cudaBindTextureToArray( imgTex, imgArray, c1 );
画素表現の指定
cudaMemcpyToArray( imgArray, 0, 0, pSrc, nBytes, cudaMemcpyHostToDevice );
転送データ量
転送元ポインタ
CPUからGPUへ転送
メモリ確保と転送（2）
テンプレートのメモリをGPU上に確保
テクスチャの定義
ハードウェア線形補間の有効化
正規化テクスチャ座標の有効化
ハードウェア線形補間の利用(1)
texture< uchar4, cudaTextureType2D, cudaReadModeNormalizedFloat > refTex;
refTex.normalized = 1;
refTex.filterMode = cudaFilterModeLinear;
ハードウェア線形補間の利用(2)
(0, 0)
(W-1, 0)
(0, H-1)
(W-1, H-1)
(0, 0)
(1, 0)
(0, 1)
(1, 1)
(0, 0)
(1, 0)
(0, 1)
(0, 0)
(1, 0)
(0, 1)
GPU上で類似度計算
__global__ void kernel( float *error, float *scale, int imgW, int imgH,
int areaW, int areaH, int maskW, int maskH, float s ) {
int i = threadIdx.x + blockDim.x * blockIdx.x; int j = threadIdx.y + blockDim.y * blockIdx.y; float err = 0.0f;
float _1_w = 1.0f / maskW; // テンプレートの幅に対するスケーリング係数 float _1_h = 1.0f / maskH; // テンプレートの高さに対するスケーリング係数 for( int n = 0 ; n < maskH ; n++ )
{
for( int m = 0 ; m < maskW ; m++ ) {
uchar4 p1 = tex2D( imgTex, i + m, j + n );
float4 p2 = tex2D( refTex, _1_w * m, _1_h * n ) * 255.0f; err += ( p1.x - p2.x ) * ( p1.x - p2.x );
err += ( p1.y - p2.y ) * ( p1.y - p2.y ); err += ( p1.z - p2.z ) * ( p1.z - p2.z ); }
}
err *= _1_w * _1_h;
if( error[ i + j * imgW ] > err ) {
error[ i + j * imgW ] = err; scale[ i + j * imgW ] = s; } }
テクスチャからデータを読み取り
SSDを計算
SSDが最小のも
のを記録
CPUによる類似度計算の制御
// スケール＆誤差を保持するGPU側のメモリ領域
thrust::device_vector< float > derror( img.size( ), 1.0e10f ); thrust::device_vector< float > dscale( img.size( ), 0 );
float s = MIN_SCALE; while( s <= MSCALE ) {
// 実行するGPUのスレッド数を指定
int W = ( int )img.width( ), H = ( int )img.height( );
int w = ( int )( ref.width( ) * s ), h = ( int )( ref.height( ) * s ); int threadNumX = 16, threadNumY = 16;
int blockNumX = ( W - w ) / threadNumX + ( ( ( W - w ) % threadNumX ) == 0 ? 0 : 1 ); int blockNumY = ( H - h ) / threadNumY + ( ( ( H - h ) % threadNumY ) == 0 ? 0 : 1 ); dim3 nThreads( threadNumX, threadNumY, 1 );
dim3 nBlocks( blockNumX, blockNumY, 1 ); // GPU側で類似度を計算する
kernel<<< nBlocks, nThreads >>>( thrust::raw_pointer_cast( &derror[ 0 ] ), thrust::raw_pointer_cast( &dscale[ 0 ] ), W, H, W - w, H - h, w, h, s ); cudaDeviceSynchronize( ); // 処理が完了するまで待機する s *= SCALE_FACTOR; // スケールを変更する } // GPUの処理結果を転送する
thrust::host_vector< float > herror = derror; thrust::host_vector< float > hscale = dscale;
GPUを用いて類似度を計算
スケールを変更して再探索
計算時間の比較（1）
CPUとGPUで計算時間を評価
使用計算機
CPU: Intel Core i7-X980（3.33GHz）
OpenMPを使用して6スレッドで並列計算
GPU: NVIDIA GeForce GTX480
マルチプロセッサ数: 15（SP数= 15×32 = 480）
Memory：1.5 GB
計算時間の比較（2）
実験パラメータ
入力画像 ：800×600画素
テンプレート：105×135画素
スケール ：0.3～1.9倍（拡大率1.2）
実験結果
CPU: 77.3 sec.
GPU: 2.2 sec.
約 35 倍の高速化
テンプレート
結果画像
まとめ（テンプレートマッチング）
テクスチャメモリを利用した高速化
キャッシュを利用した効率的なメモリアクセス
解像度に依存しないメモリアクセス
ハードウェア線形補間
正規化テクスチャ座標
CPUとGPUで計算時間を比較
CPU: 77.3 sec.
GPU: 2.2 sec.
約 35 倍の高速化を実現
類似度計算の打ち切り処理を導入
各スケールで計算した類似度の最大値をGPUで計算
Parallel Reduction アルゴリズムを利用
スケールの変更
画
像
読
み
込
み
CPU
メ
モ
リ
確
保
GPU
類
似
度
計
算
GPU
ス
レ
ッ
ド
同
期
CPU
デ
ー
タ
転
送
CPU GPU
デ
ー
タ
転
送
GPU CPU
最
大
類
似
度
計
算
GPU
処理の流れ
GPU上で類似度計算
__global__ void kernel( float *error, float *scale, int imgW, int imgH,
int areaW, int areaH, int maskW, int maskH, float s, float maxerr ) {
int i = threadIdx.x + blockDim.x * blockIdx.x; int j = threadIdx.y + blockDim.y * blockIdx.y; float err = 0.0f;
float _1_w = 1.0f / maskW; // テンプレートの幅に対するスケーリング係数 float _1_h = 1.0f / maskH; // テンプレートの高さに対するスケーリング係数 for( int n = 0 ; n < maskH ; n++ )
{
for( int m = 0 ; m < maskW ; m++ ) {
uchar4 p1 = tex2D( imgTex, i + m, j + n );
float4 p2 = tex2D( refTex, _1_w * m, _1_h * n ) * 255.0f; err += ( p1.x - p2.x ) * ( p1.x - p2.x );
err += ( p1.y - p2.y ) * ( p1.y - p2.y ); err += ( p1.z - p2.z ) * ( p1.z - p2.z ); }
if( maxerr < err *_1_w * _1_h ) { break; } } ... 以下同じ ... }
前スケールの探索時における
SSDが最小値より大きい場合
は計算を打ち切り
Parallel Reductionによる最小値探索
GPU上のメモリから最小値を計算
GPUのスレッドが並列に値の比較処理を実行
2
GPU上のデータ
同期
11 4 27 25 5 13 6 9 20 12 14 7 2 19 3 15
4
25
5
6
12
7
2
3
4
5
7
2
4
2
2
同期
同期
Thrustライブラリの reduce アルゴリズムを利用
thrust::reduce( 先頭, 末尾, 初期値, thrust::minimum<T>( ) );
Parallel Reductionによる最小値探索
GPU上のメモリから最小値を計算
GPUのスレッドが並列に値の比較処理を実行
2
GPU上のデータ
同期
11 4 27 25 5 13 6 9 20 12 14 7 2 19 3 15
4
25
5
6
12
7
2
3
4
5
7
2
4
2
2
同期
同期
先頭を指すイテレータ
末尾を指すイテレータ
CPUによる類似度計算の制御
// スケール＆誤差を保持するGPU側のメモリ領域
thrust::device_vector< float > derror( img.size( ), 1.0e10f ); thrust::device_vector< float > dscale( img.size( ), 0 );
float s = MIN_SCALE, maxerr = 1.0e10f; while( s <= MSCALE )
{
// 実行するGPUのスレッド数を指定
int W = ( int )img.width( ), H = ( int )img.height( );
int w = ( int )( ref.width( ) * s ), h = ( int )( ref.height( ) * s ); int threadNumX = 16, threadNumY = 16;
int blockNumX = ( W - w ) / threadNumX + ( ( ( W - w ) % threadNumX ) == 0 ? 0 : 1 ); int blockNumY = ( H - h ) / threadNumY + ( ( ( H - h ) % threadNumY ) == 0 ? 0 : 1 ); dim3 nThreads( threadNumX, threadNumY, 1 );
dim3 nBlocks( blockNumX, blockNumY, 1 ); // GPU側で類似度を計算する
kernel<<< nBlocks, nThreads >>>( thrust::raw_pointer_cast( &derror[ 0 ] ), thrust::raw_pointer_cast( &dscale[ 0 ] ),
W, H, W - w, H - h, w, h, s, maxerr );
cudaDeviceSynchronize( ); // 処理が完了するまで待機する // 誤差の最大値を取得する
maxerr = thrust::reduce( derror.begin( ), derror.end( ), 1.0e10f, thrust::minimum< float >( ) ); s *= SCALE_FACTOR; // スケールを変更する } // GPUの処理結果を転送する ...
