GPUによる画像処理

授業項目
並列計算の概念（プロセスとスレッド）
9. GPUのアーキテクチャ
10. GPGPUのプログラム構造
11. GPUでの並列プログラミング（ベクトル和）
12. GPUによる画像処理
13. GPUのメモリ階層
14. GPGPU組込開発環境
15. GPGPU開発環境（OpenCL）
8.
506
GPGPU実践基礎工学
2015/11/25
第12回 GPUによる画像処理
長岡技術科学大学
電気電子情報工学専攻
出川智啓
今回の内容

画像処理

ネガティブ処理，画像反転，空間フィルタ，モザイク処理

画像処理を簡略化した1次元的な処理

2次元的な並列処理



2次元配列の表現
ベクトル和
モノクロ画像処理
508
GPGPU実践基礎工学
2015/11/25
画像処理

デジタル画像を別の画像に変換

処理の種類

画像の拡大・縮小

画像の反転

フィルタ

509
2値化，ぼかしなど
GPGPU実践基礎工学
2015/11/25
画像処理を行うプログラムの流れ

読み込む画像を準備

対象ファイルをオープン

画像の読込



画像ファイルの種類（bmp,jpg,png等）に応じて読み込む方
法を選択
画像ファイルからヘッダ情報を読み込み，様々な情報を取得
メモリを確保し，画素情報を読込

定められたアルゴリズムで画像を処理

処理された画像の出力
510
GPGPU実践基礎工学
2015/11/25
画像処理

画像を読むプログラムを作るだけで一苦労

画像処理が主目的ではなく，GPUの使い方を学習するこ
とが目的

画像をモノクロに限定




プログラム中で画像に相当するデータを生成
ターミナルにテキスト出力（横の位置，縦の位置，色情報）
画像の表示にはgnuplotを利用
モノクロ画像に対して簡単な処理を実行
511
GPGPU実践基礎工学
2015/11/25
取り扱うモノクロ画像

型はunsigned char


1バイトの整数型
取り扱える範囲は0~255

0を黒，255を白とみなす

1次元配列で取り扱う


512
画像は2次元なので，2次元配列で取り扱う方が自然
mallocやcudaMallocでは2次元配列を確保できない
GPGPU実践基礎工学
2015/11/25
取り扱うモノクロ画像の内容
x
unsigned char
picture[] = { 255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,
255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,
255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,
255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,
255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,
255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,
255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,
255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,
y 255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,
255,255,255,255,255,255,255,255,255, 0,255,255,255,255,255,255,
255,255,255,255,255,255,255,255,255, 0, 0,255,255,255,255,255,
255,255,255,255,255,255,255,255,255, 0, 0, 0,255,255,255,255,
255,255,255,255,255,255,255,255,255, 0, 0, 0, 0,255,255,255,
255,255,255,255,255,255,255,255,255, 0, 0, 0, 0, 0,255,255,
255,255,255,255,255,255,255,255,255, 0, 0, 0, 0, 0, 0,255,
255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255};
513
GPGPU実践基礎工学
2015/11/25
取り扱うモノクロ画像の内容

gnuplotで表示


左下を原点として表示
画像処理では左上が原点
y
x
514
GPGPU実践基礎工学
2015/11/25
画像処理

ネガティブ処理

反転


空間フィルタ


水平反転，垂直反転
ぼかし，輪郭抽出
モザイク処理
515
GPGPU実践基礎工学
2015/11/25
画像処理の簡略化（1次元化）

いきなり画像処理は難しい

1次元の配列に対して，画像処理と同じような処理を実行

どのような処理を行うかを確認
516

ネガティブ処理

水平反転

空間フィルタ（ぼかし，輪郭抽出）

モザイク処理
GPGPU実践基礎工学
原画像もどき



幅 WIDTH=64 (0~63)
[3*WIDTH/16,9*WIDTH/16)の範囲を白(255)，それ
以外を黒(0)
変数名


原画像pのピクセル値
画素値

原画像 p[]
処理画像 filtered[]
横位置
517
GPGPU実践基礎工学
2015/11/25
ネガティブ処理

色の反転処理

白(255)を黒(0)に
原画像p[]
黒(0)を白(255)に
処理画像
filtered[]
画素値
画素値

黒(0)を白(255)，白(255)を黒(0)に反転
(白−各画素の値)を新しい画素の値とする
横位置
518
GPGPU実践基礎工学
2015/11/25
ネガティブ処理

色の反転処理


黒(0)を白(255)，白(255)を黒(0)に反転
(白−各画素の値)を新しい画素の値とする
void negative(unsigned char *p, unsigned char *filtered){
int i;
for (i = 0; i<WIDTH; i++){
filtered[i] = WHITE ‐ p[i];
}
}
519
GPGPU実践基礎工学
2015/11/25
水平反転

画像の幾何的位置を反転

画素値
画素値

水平反転 i→(WIDTH‐1)‐i
p[i]→filtered[(WIDTH‐1)‐i]
横位置1の画素値を62番目の画素値に
横位置35の画素値を28番目の画素値に
原画像p[]
処理画像
filtered[]
横位置
520
GPGPU実践基礎工学
2015/11/25
水平反転

画像の幾何的位置を反転


水平反転 i→(WIDTH‐1)‐i
p[i]→filtered[(WIDTH‐1)‐i]
void hreflect(unsigned char *p, unsigned char *filtered){
int i;
for (i = 0; i<WIDTH; i++){
//(WIDTH‐1)‐iは反転後の横位置
filtered[(WIDTH ‐ 1) ‐ i] = p[i];
}
}
521
GPGPU実践基礎工学
2015/11/25
空間フィルタ

ある画素とその周囲の画素を使って処理



処理の仕方を規定したフィルタカーネルを定義
フィルタカーネルは1次元配列で表現，引数として関数に渡す
ぼかし（平均）フィルタ
1/3

1/3
1/3
float blur[3] ={1.0f/3.0f,1.0f/3.0f,1.0f/3.0f};
輪郭抽出（ラプラシアン）フィルタ
1
522
‐2
1
float laplacian[9] ={ 1.0f,‐2.0f, 1.0f};
GPGPU実践基礎工学
2015/11/25
空間フィルタ

ある画素とその周囲の画素を使って処理
画素値

ある位置iの左右(i‐1,i+1)の計3画素の値を読み，それぞ
れフィルタカーネルの各値との積を計算し，合計を取った値を
処理画像のi番目の画素とする
原画像p[]
ぼかしフィルタカーネル
1/3
1/3
1/3
255/3+255/3+255/3=255
画素値
処理画像
filtered[]
横位置
523
GPGPU実践基礎工学
2015/11/25
空間フィルタ

ある画素とその周囲の画素を使って処理
画素値

ある位置iの左右(i‐1,i+1)の計3画素の値を読み，それぞ
れフィルタカーネルの各値との積を計算し，合計を取った値を
処理画像のi番目の画素とする
原画像p[]
ぼかしフィルタカーネル
1/3
1/3
1/3
255/3+255/3+255/3=255
画素値
処理画像
filtered[]
横位置
524
GPGPU実践基礎工学
2015/11/25
空間フィルタ

ある画素とその周囲の画素を使って処理
画素値

ある位置iの左右(i‐1,i+1)の計3画素の値を読み，それぞ
れフィルタカーネルの各値との積を計算し，合計を取った値を
処理画像のi番目の画素とする
原画像p[]
ぼかしフィルタカーネル
1/3
1/3
1/3
255/3+255/3+255/3=255
画素値
処理画像
filtered[]
横位置
525
GPGPU実践基礎工学
2015/11/25
空間フィルタ

ある画素とその周囲の画素を使って処理
ある位置iの左右(i‐1,i+1)の計3画素の値を読み，それぞ
れフィルタカーネルの各値との積を計算し，合計を取った値を
処理画像のi番目の画素とする
画素値

原画像p[]
ぼかしフィルタカーネル
1/3
1/3
1/3
0/3+255/3+255/3=170
画素値
処理画像
filtered[]
横位置
526
GPGPU実践基礎工学
2015/11/25
空間フィルタ

ある画素とその周囲の画素を使って処理

forループでiの値を変化させ，p[i‐1],p[i],p[i+1]の値
から処理後の画素値filtered[i]を計算
void boxfilter
(unsigned char *p, unsigned char *filtered, float *filter){
int i, result=BLACK;
for (i = 1; i<WIDTH ‐ 1; i++){//端の画素は処理をしない
result = filter[0]*p[i‐1]+filter[1]*p[i]+filter[2]*p[i+1];
//フィルタ後の値が負になれば0に切り上げ，255を超えたら255に収める
if (result<BLACK) result = 0;
if (result>WHITE) result = WHITE;
filtered[i] = (unsigned char)result;
}
}
527
GPGPU実践基礎工学
2015/11/25
空間フィルタ

輪郭抽出フィルタの結果

画素値
画素値

白と黒の境界でのみ値が255になり，それ以外では0
両端はフィルタ処理できないため，初期値がそのまま使われる
原画像p[]
白と黒の境界でのみ
画素値が0より大きく
なる
処理画像
filtered[]
横位置
528
GPGPU実践基礎工学
2015/11/25
モザイク処理

画像を小領域に分け，その領域を全て同じ色に置換

画素値
原画像p[]
画素値

小領域内の全画素を，小領域内の画素の平均値に置換
小領域を移動するループ内に，小領域内の画素の平均値を
計算するループを設ける
処理画像
filtered[]
小領域
横位置
529
GPGPU実践基礎工学
2015/11/25
モザイク処理

画像を小領域に分け，その領域を全て同じ色に置換

画素値
原画像p[]
画素値

小領域内の全画素を，小領域内の画素の平均値に置換
小領域を移動するループ内に，小領域内の画素の平均値を
計算するループを設ける
処理画像
filtered[]
横位置
530
GPGPU実践基礎工学
2015/11/25
モザイク処理

画像を小領域に分け，その領域を全て同じ色に置換

画素値
原画像p[]
画素値

小領域内の全画素を，小領域内の画素の平均値に置換
小領域を移動するループ内に，小領域内の画素の平均値を
計算するループを設ける
処理画像
filtered[]
小領域
横位置
531
GPGPU実践基礎工学
2015/11/25
モザイク処理

画像を小領域に分け，その領域を全て同じ色に置換

画素値
原画像p[]
画素値

小領域内の全画素を，小領域内の画素の平均値に置換
小領域を移動するループ内に，小領域内の画素の平均値を
計算するループを設ける
処理画像
filtered[]
横位置
532
GPGPU実践基礎工学
2015/11/25
モザイク処理

画像を小さな領域に分け，その領域を全て同じ色に置換

小領域内で画素の平均値を計算し，処理画像の小領域内の
全画素をその平均値とする
void mosaic(unsigned char *p,unsigned char *filtered,int mosaicSize){
int i,isub,average;
for (i = 0; i<WIDTH; i += mosaicSize){//小領域を移動するループ
//小領域内の画素の平均値を計算
average = 0;
for(isub = 0; isub<mosaicSize; isub++)
average += p[(i+isub)]; average /= mosaicSize; //小領域内の全画素を平均値で塗りつぶす
for(isub = 0; isub<mosaicSize; isub++)
filtered[(i+isub)] = (unsigned char)average;
}
}
533
GPGPU実践基礎工学
2015/11/25
処理の流れ

変数を宣言し，mallocでメモリを確保



原画像pを作成


pをfilteredにコピーしておく
処理を実行


原画像を格納する変数
p
処理後の画像を格納する変数 filtered
関数を作成し，pとfilteredを渡す
filteredの内容を画面に出力
534
GPGPU実践基礎工学
2015/11/25
メイン関数（CPU版）
#include<stdio.h>
#include<stdlib.h>
#define WHITE (255)
#define BLACK (0)
#define WIDTH (64)
#define Nbytes (WIDTH*sizeof(unsigned char))
void create(unsigned char *);
void copy(unsigned char *,unsigned char *);
void print(unsigned char *);
void negative(unsigned char *, unsigned char *);
void hreflect(unsigned char *, unsigned char *);
void boxfilter(unsigned char *, unsigned char *, float *);
void mosaic(unsigned char *, unsigned char *, int);
int main(void){
unsigned char *p= (unsigned char *)malloc(Nbytes);
unsigned char *filtered= (unsigned char *)malloc(Nbytes);
//ここで空間フィルタのカーネルを宣言
535
imageproc1d.c
GPGPU実践基礎工学
2015/11/25
メイン関数（CPU版）
create(p);
copy(p,filtered);
//ここで処理を行い，結果をfilteredに格納
//画面に各画素の値を表示
print(filtered);
return 0;
}
//画像の内容をコピー
void copy(unsigned char *src, unsigned char *dst){
int i;
for (i = 0; i<WIDTH; i++)
dst[i] = src[i];
}
//画像の内容を画面に出力
void print(unsigned char *p){
int i;
for (i = 0; i<WIDTH; i++)
printf("%d %d¥n", i, p[i]);
}
536
GPGPU実践基礎工学
imageproc1d.c
2015/11/25
画像の作成

[3*WIDTH/16,9*WIDTH/16)の範囲を白(255)，それ
以外を黒(0)
void create(unsigned char *p){
int i, x_origin;
for (i = 0; i<WIDTH; i++)
p[i] = WHITE;
x_origin = 3 * WIDTH / 16;
for (i = 0; i<6 * WIDTH / 16; i++)
p[i + x_origin] = BLACK;
画素値
}
3WIDTH/16
537
6WIDTH/16
GPGPU実践基礎工学
横位置
2015/11/25
結果の画面表示とファイル保存

print関数が画面に各画素の座標と値を表示
0 255
1 255
:
62 255
63 255

画面出力のファイル表示（出力のリダイレクトを利用）


>を利用して保存，結果は表計算等で確認
$ コマンド > ファイル名
$ ./a.out > pic.txt
$ cat pic.txt
0 255
1 255
:
62 255
63 255
538
GPGPU実践基礎工学
2015/11/25
GPUへの移植の方針

現画像の作成はCPUで行い，GPUへ転送

1スレッドが1画素の処理を実行


空間フィルタ


forループを排除し，配列添字とスレッド番号を対応
端の画素の処理ができないので今回は処理しない
モザイク


領域の大きさは1ブロックと同じとする
領域内の画素の平均値の計算は非常に難易度が高い

539
各ブロックにおけるthreadIdx.x=0のスレッドが平均を計算して処
理画像に書き込む
GPGPU実践基礎工学
2015/11/25
メイン関数（GPU版）
#include<stdio.h>
#include<stdlib.h>
#define WHITE (255)
#define BLACK (0)
#define WIDTH (64)
#define Nbytes (WIDTH*sizeof(unsigned char))
#define THREAD 16
#define BLOCK (WIDTH/THREAD)
void create(unsigned char *);
void print(unsigned char *);
int main(void){
unsigned char *p= (unsigned char *)malloc(Nbytes);
unsigned char *dev_p, *dev_filtered;
cudaMalloc( (void **)&dev_p, Nbytes);
cudaMalloc( (void **)&dev_filtered, Nbytes);
imageproc1d.cu
540
GPGPU実践基礎工学
2015/11/25
メイン関数（GPU版）
//ここで空間フィルタのカーネルを宣言
float *filter;
cudaMalloc( (void **)&filter, sizeof(float)*3);
cudaMemcpy(filter, フィルタのカーネル, sizeof(float)*3, cudaMemcpyHostToDevice);
create(p); //CPUで画像を生成してGPUへ送った後，変数dev_filteredにコピーしておく
cudaMemcpy(dev_p, p, Nbytes, cudaMemcpyHostToDevice);
cudaMemcpy(dev_filtered, dev_p, Nbytes, cudaMemcpyDeviceToDevice);
//ここで処理を行い，結果をdev_filteredに格納
//dev_filteredの内容をCPUへ送る
unsigned char *filtered = (unsigned char *)malloc(Nbytes);
cudaMemcpy(filtered, dev_filtered, Nbytes, cudaMemcpyDeviceToHost);
free(p);
free(filtered);
cudaFree(dev_p);
cudaFree(dev_filtered);
cudaFree(filter);
imageproc1d.cu
541
GPGPU実践基礎工学
2015/11/25
メイン関数（GPU版）
return 0;
}
//原画像の生成
void create(unsigned char *p){
int i, x_origin;
for (i = 0; i<WIDTH; i++)
p[i] = WHITE;
x_origin = 3 * WIDTH / 16;
for (i = 0; i<6 * WIDTH / 16; i++)
p[i + x_origin] = BLACK;
}
//画像の内容を画面に出力
void print(unsigned char *p){
int i;
for (i = 0; i<WIDTH; i++)
printf("%d %d¥n", i, p[i]);
}
imageproc1d.cu
542
GPGPU実践基礎工学
2015/11/25
ネガティブ処理

1スレッドが1画素の値を計算

ベクトル和のようにスレッド番号と配列添字の対応を計算
画素値
原画像p[]
画素値
スレッド番号 0 1 2 3 ･･･
処理画像
filtered[]
横位置
543
GPGPU実践基礎工学
2015/11/25
ネガティブ処理
void negative(unsigned char *p, unsigned char *filtered){
int i;
for (i = 0; i<WIDTH; i++)
filtered[i] = WHITE ‐ p[i];
}
imageproc1d.c
544
GPGPU実践基礎工学
2015/11/25
ネガティブ処理
__global__ void negative(unsigned char *p, unsigned char *filtered){
int i;
i = blockIdx.x*blockDim.x + threadIdx.x; //スレッド数と配列添字の対応
filtered[i] = WHITE ‐ p[i];
}
imageproc1d.cu
545
GPGPU実践基礎工学
2015/11/25
水平反転

1スレッドが1画素の値を反転


ベクトル和のようにスレッド番号と配列添字の対応を計算
書込先が変わる
画素値
原画像p[]
画素値
スレッド番号 0 1 2 3 ･･･
処理画像
filtered[]
横位置35の画素値を28番目の画素値に
横位置
546
GPGPU実践基礎工学
2015/11/25
水平反転
void hreflect(unsigned char *p, unsigned char *filtered){
int i;
for (i = 0; i<WIDTH; i++)
//(WIDTH‐1)‐iは反転後の横位置
filtered[(WIDTH ‐ 1) ‐ i] = p[i];
}
imageproc1d.c
547
GPGPU実践基礎工学
2015/11/25
水平反転
__global__ void hreflect(unsigned char *p, unsigned char *filtered){
int i;
i = blockIdx.x*blockDim.x + threadIdx.x;
//(WIDTH‐1)‐iは反転後の横位置
filtered[(WIDTH‐1)‐i] = p[i];
}
imageproc1d.cu
548
GPGPU実践基礎工学
2015/11/25
空間フィルタ

1スレッドが1画素の値を計算


ある位置iの画素だけではなく左右(i‐1,i+1)の画素も参照
全スレッドが同じフィルタカーネルを参照
原画像p[]
1/3
1/3
1/3
ぼかしフィルタカーネル
処理画像
filtered[]
画素値
画素値
スレッド番号 0 1 2 3 ･･･
横位置
549
GPGPU実践基礎工学
2015/11/25
空間フィルタ
void boxfilter
(unsigned char *p, unsigned char *filtered, float *filter){
int i, result=BLACK;
for (i = 1; i<WIDTH ‐ 1; i++){//端の画素は処理をしない
result = filter[0]*p[i‐1]+filter[1]*p[i]+filter[2]*p[i+1];
//フィルタ後の値が負になれば0に切り上げ，255を超えたら255に収める
if (result<BLACK) result = 0;
if (result>WHITE) result = WHITE;
filtered[i] = (unsigned char)result;
}
}
imageproc1d.c
550
GPGPU実践基礎工学
2015/11/25
空間フィルタ
__global__ void boxfilter
(unsigned char *p, unsigned char *filtered, float *filter){
int i, result = BLACK;
i = blockIdx.x*blockDim.x + threadIdx.x;
if(0<i && i<WIDTH‐1) //端の画素は処理をしないため，ifで処理を分岐
result = filter[0]*p[i‐1]+filter[1]*p[i]+filter[2]*p[i+1];
//フィルタ後の値が負になれば0に切り上げ，255を超えたら255に収める
if(result<BLACK) result = 0;
if(result>WHITE) result = WHITE;
filtered[i] = (unsigned char)result;
}
imageproc1d.cu
551
GPGPU実践基礎工学
2015/11/25
モザイク処理

1ブロック内の1スレッドのみが処理


1スレッドがブロック内の全画素の平均を計算
平均の画素値を処理画像に書込
画素値
原画像p[]
画素値
スレッド番号 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7 ･･･
処理画像
filtered[]
横位置
552
GPGPU実践基礎工学
2015/11/25
モザイク処理
void mosaic(unsigned char *p,unsigned char *filtered,int mosaicSize){
int i,isub,average;
for (i = 0; i<WIDTH; i += mosaicSize){//小領域を移動するループ
//小領域内の画素の平均値を計算
average = 0;
for(isub = 0; isub<mosaicSize; isub++)
average += p[(i+isub)]; average /= mosaicSize; //小領域内の全画素を平均値で塗りつぶす
for(isub = 0; isub<mosaicSize; isub++)
filtered[(i+isub)] = (unsigned char)average;
}
}
imageproc1d.c
553
GPGPU実践基礎工学
2015/11/25
モザイク処理
__global__ void mosaic
(unsigned char *p, unsigned char *filtered, int mosaicSize){
int i,isub,average;
i = blockIdx.x*blockDim.x + threadIdx.x;
if(threadIdx.x == 0){//ブロック内の1スレッドのみが処理
//小領域内の画素の平均値を計算
average = 0;
for(isub = 0; isub<mosaicSize; isub++)
average += p[(i+isub)];
average /= mosaicSize;
//小領域内の全画素を平均値で塗りつぶす
for(isub = 0; isub<mosaicSize; isub++)
filtered[(i+isub)] = (unsigned char)average;
}
}
imageproc1d.cu
554
GPGPU実践基礎工学
2015/11/25
ネガティブ処理

色の反転処理

黒を白，白を黒に反転
原画像
555
処理画像
GPGPU実践基礎工学
2015/11/25
画像反転

画像の幾何的位置を反転

画像（画素）の配列添字を変えて出力
原画像
556
処理画像
GPGPU実践基礎工学
2015/11/25
空間フィルタ

ある画素とその周囲の画素を使って処理


処理の仕方を規定したカーネルを定義
カーネルを変更することで様々な処理が可能
原画像
557
輪郭抽出
GPGPU実践基礎工学
2015/11/25
フィルタ処理

2次元的な処理


あるピクセルの周囲情報を利用して効果を計算
画像を矩形領域で分割した2次元的な並列処理が自然
原画像
558
輪郭抽出
フィルタ
（カーネル）
a
b
c
d
e
f
g
h
i
0
1
0
1
‐4
1
0
1
0
GPGPU実践基礎工学
= b+d‐4e+f+h
2015/11/25
モザイク処理

画像を小さな領域に分け，その領域を全て同じ色に置換

領域内のある1画素の値か，領域内の画素の平均値にするこ
とが多い
原画像
559
処理画像
GPGPU実践基礎工学
2015/11/25
2次元配列

変数名の後ろに[]を2個付けて宣言
unsigned char picture[Nx][Ny];

1次元配列はベクトル，2次元配列は行列や画像
i
i
matrix[Nx][Ny]
vector[Nx]
Ny=8
Nx=8
j
Nx=8
560
GPGPU実践基礎工学
2015/11/25
2次元配列の1次元配列表現

mallocやcudaMallocは1次元配列を宣言



2次元配列を宣言できない
1次元配列を宣言し，2次元配列的に読み書き
アドレス空間は1次元


2次元配列でも1次元のアドレスで表現
1次元目が連続か，2次元目が連続かが言語によって異なる
1次元目が連続
2次元目が連続
i
i
561
matrix[Nx][Ny]
matrix[Nx][Ny]
j
j
GPGPU実践基礎工学
2015/11/25
2次元配列の1次元配列表現

C言語は2次元目が連続
picture[Nx][Ny]
i
2次元配列の場合
picutre[1][2] = 11

1次元配列の場合
picture[_______] = 11
1*8 + 2
9
10
11
12
Ny=8

1
2
3
4
5
6
7
j
8
Nx=8
picture[i][j]  picture[i*Ny+j]
562
GPGPU実践基礎工学
2015/11/25
ベクトル和C=A+B（2次元版）

1次元のベクトル和とほぼ同じ


並列化が容易
for文の書く順番で実行速度が変化
c[][]
563
a[][]
GPGPU実践基礎工学
b[][]
2015/11/25
ベクトル和（2次元版）
#define Nx (1024)
#define Ny (1024)
//iとjのループを入れ替えるとどうなるか？
for(i=0;i<Nx;i++){
for(j=0;j<Nx;j++){
c[i][j] = a[i][j]+b[i][j];
}
}
void init(float a[Nx][Ny], float b[Nx][Ny],
float c[Nx][Ny]){
}
int i,j;
//iとjのループを入れ替えるとどうなるか？
for(i=0;i<Nx;i++){
for(j=0;j<Nx;j++){
a[i][j] = 1.0;
b[i][j] = 2.0;
c[i][j] = 0.0;
}
}
int main(void){
float a[Nx][Ny],b[Nx][Ny],c[Nx][Ny];
init(a, b, c);
add(a, b, c);
return 0;
}
}
void add(float a[Nx][Ny], float b[Nx][Ny],
float c[Nx][Ny]){
int i,j;
564
vectoradd2d.c
GPGPU実践基礎工学
2015/11/25
ベクトル和（2次元→1次元版）
for(j=0;j<Ny;j++){
ij = i*Ny+j;
//ij = i+Nx*j;にするとどうなるか？
c[ij] = a[ij]+b[ij];
}
#include<stdlib.h>
#define Nx (1024)
#define Ny (1024)
#define Nbytes (Nx*Ny*sizeof(float))
}
void init(float *a, float *b, float *c){
}
int i,j,ij;
for(i=0;i<Nx;i++){
for(j=0;j<Ny;j++){
ij = i*Ny+j;
//ij = i+Nx*j;にするとどうなるか？
a[ij] = 1.0;
b[ij] = 2.0;
c[ij] = 0.0;
}
}
int main(void){
float *a,*b,*c;
a = (float *)malloc(Nbytes);
b = (float *)malloc(Nbytes);
c = (float *)malloc(Nbytes);
init(a, b, c);
add(a, b, c);
}
void add(float *a, float *b, float *c){
int i,j,ij;
for(i=0;i<Nx;i++){
565
free(a);
free(b);
free(c);
return 0;
}
GPGPU実践基礎工学
vectoradd2d_1darray.c
2015/11/25
GPUによる2次元的な並列処理

ブロックとスレッドを2次元的に設定


566
1スレッドが配列の1要素（画像の1ピクセル）を処理
スレッドを2次元的に配置してブロックを構成
GPGPU実践基礎工学
2015/11/25
GPUの並列化の階層


グリッド－ブロック－スレッドの3階層
各階層の情報を参照できる変数


グリッド（Grid）


gridDim
グリッド内にあるブロックの数
ブロック（Block）



x,y,zをメンバにもつdim3型構造体
blockIdx
blockDim
ブロックに割り当てられた番号
ブロック内にあるスレッドの数
スレッド（Thread）

567
threadIdx
スレッドに割り当てられた番号
GPGPU実践基礎工学
2015/11/25
Hello Threads（2次元版）

<<< , >>>の中にどのように数字を書くか
#include<stdio.h>
__global__ void hello(){
printf("gridDim.x=%d,blockIdx.x=%d,blockDim.x=%d,threadIdx.x=%d¥n",
gridDim.x, blockIdx.x, blockDim.x, threadIdx.x);
printf("gridDim.y=%d,blockIdx.y=%d,blockDim.y=%d,threadIdx.y=%d¥n",
gridDim.y, blockIdx.y, blockDim.y, threadIdx.y);
}
int main(void){
hello<<<?,?>>>();
cudaDeviceSynchronize();
return 0;
}
hellothreads2d.cu
568
GPGPU実践基礎工学
2015/11/25
Hello Threads（2次元版）

1次元と同様<<<2,4>>>等として実行

実行結果（画面出力）
gridDim.x=2, blockIdx.x=0, blockDim.x=4, threadIdx.x=0
gridDim.x=2, blockIdx.x=0, blockDim.x=4, threadIdx.x=1
gridDim.x=2, blockIdx.x=0, blockDim.x=4, threadIdx.x=2
gridDim.x=2, blockIdx.x=0, blockDim.x=4, threadIdx.x=3
gridDim.x=2, blockIdx.x=1, blockDim.x=4, threadIdx.x=0
gridDim.x=2, blockIdx.x=1, blockDim.x=4, threadIdx.x=1
gridDim.x=2, blockIdx.x=1, blockDim.x=4, threadIdx.x=2
gridDim.x=2, blockIdx.x=1, blockDim.x=4, threadIdx.x=3
gridDim.y=1, blockIdx.y=0, blockDim.y=1, threadIdx.y=0
gridDim.y=1, blockIdx.y=0, blockDim.y=1, threadIdx.y=0
gridDim.y=1, blockIdx.y=0, blockDim.y=1, threadIdx.y=0
gridDim.y=1, blockIdx.y=0, blockDim.y=1, threadIdx.y=0
gridDim.y=1, blockIdx.y=0, blockDim.y=1, threadIdx.y=0
gridDim.y=1, blockIdx.y=0, blockDim.y=1, threadIdx.y=0
gridDim.y=1, blockIdx.y=0, blockDim.y=1, threadIdx.y=0
gridDim.y=1, blockIdx.y=0, blockDim.y=1, threadIdx.y=0
569
GPGPU実践基礎工学
並列度の指定が
できていない
2015/11/25
2次元的な並列度の指定

<<<,>>>の中にどのように数字を書くか



1次元の場合は数字を書くことができた
2次元，3次元は数字を並べて書くことができない
dim3型構造体を利用
int main(void){
dim3 block(2,4,1), thread(4,2,1); ・・・ dim3型変数block, hello<<<block, thread>>>();
threadを利用
hello<<<dim3(2,4,1), dim3(4,2,1)>>>(...);
・・・あるいは直接dim3型として記述
}
570
GPGPU実践基礎工学
2015/11/25
Hello Threads（2次元版）

実行結果（画面出力）
gridDim.x=2, blockIdx.x=0, blockDim.x=4, threadIdx.x=0
gridDim.x=2, blockIdx.x=0, blockDim.x=4, threadIdx.x=1
gridDim.x=2, blockIdx.x=0, blockDim.x=4, threadIdx.x=2
gridDim.x=2, blockIdx.x=0, blockDim.x=4, threadIdx.x=3
gridDim.x=2, blockIdx.x=0, blockDim.x=4, threadIdx.x=0
gridDim.x=2, blockIdx.x=0, blockDim.x=4, threadIdx.x=1
gridDim.x=2, blockIdx.x=0, blockDim.x=4, threadIdx.x=2
gridDim.x=2, blockIdx.x=0, blockDim.x=4, threadIdx.x=3
:(中略)
gridDim.y=4, blockIdx.y=0, blockDim.y=2, threadIdx.y=0
gridDim.y=4, blockIdx.y=0, blockDim.y=2, threadIdx.y=0
gridDim.y=4, blockIdx.y=0, blockDim.y=2, threadIdx.y=0
gridDim.y=4, blockIdx.y=0, blockDim.y=2, threadIdx.y=0
gridDim.y=4, blockIdx.y=0, blockDim.y=2, threadIdx.y=1
gridDim.y=4, blockIdx.y=0, blockDim.y=2, threadIdx.y=1
gridDim.y=4, blockIdx.y=0, blockDim.y=2, threadIdx.y=1
gridDim.y=4, blockIdx.y=0, blockDim.y=2, threadIdx.y=1
:
571
GPGPU実践基礎工学
対応
対応
2015/11/25
各スレッドが異なるi,jを参照するには

CUDAでは並列化に階層がある


全体の領域（グリッド）をブロックに分割
ブロックの中をスレッドに分割
<<<dim3(2,4,1), dim3(4,2,1)>>>
ブロックの数
1ブロックあたりの
スレッドの数
x方向ブロック数×1ブロックあたりのスレッドの数
×y方向ブロック数×1ブロックあたりのスレッドの数＝総スレッド数
総スレッド数=2×4×4×2=64
572
GPGPU実践基礎工学
2015/11/25
各スレッドが異なるi,jを参照するには

Nx=8, Ny=8, x,y方向スレッド数4，ブロック数2


gridDim.x=2, gridDim.y=2
blockDim.x=4,blockDim.y=4
threadIdx.x,threadIdx.y
2 3
4
(0,0)(1,0)(2,0)(3,0)(0,0)
5 6
7
・・・
(0,1)(1,1)(2,1)(3,1)
(0,3)(1,3)(2,3)(3,3)
・・・
(0,0)
・・・
5 6
・・・
・・・
GPGPU実践基礎工学
(3,3)
7
(3,3)
573
4
(0,0)
(3,3)
2 3
・・・
(0,2)(1,2)(2,2)(3,2)
j= 0 1
i= 0 1
2015/11/25
各スレッドが異なるi,jを参照するには

Nx=8, Ny=8, x,y方向スレッド数4，ブロック数2
gridDim.x=2
574
blockIdx.y=1
gridDim.y=2
threadIdx.x,threadIdx.y
(0,0)(1,0)(2,0)(3,0)(0,0)
(0,1)(1,1)(2,1)(3,1)
(0,2)(1,2)(2,2)(3,2)
(0,3)(1,3)(2,3)(3,3)
(0,0)
(3,3)
blockDim.y=4
blockIdx.y=0
blockIdx.x=0 blockIdx.x=1
(0,0)
(3,3)
(3,3)
blockDim.x=4
GPGPU実践基礎工学
2015/11/25
各スレッドが異なるi,jを参照するには

Nx=8, Ny=8, x,y方向スレッド数4，ブロック数2


i = blockIdx.x*blockDim.x + threadIdx.x
j = blockIdx.y*blockDim.y + threadIdx.y
threadIdx.x,threadIdx.y
block(1,0)
4 5 6 7
(0,0)(1,0)(2,0)(3,0)(0,0)
(0,1)(1,1)(2,1)(3,1)
(0,3)(1,3)(2,3)(3,3)
4
(0,0)
(3,3)
2 3
(0,2)(1,2)(2,2)(3,2)
j= 0 1
block(0,0)
i= 0 1 2 3
(0,0)
5 6
block(0,1)
575
(3,3)
7
(3,3)
block(1,1)
GPGPU実践基礎工学
2015/11/25
ベクトル和（1次元配列版）
ij = i*Ny+j;
c[ij] = a[ij]+b[ij];
#include<stdlib.h>
#define Nx (1024)
#define Ny (1024)
#define Nbytes (Nx*Ny*sizeof(float))
}
}
}
void init(float *a, float *b, float *c){
int main(void){
float *a,*b,*c;
a = (float *)malloc(Nbytes);
b = (float *)malloc(Nbytes);
c = (float *)malloc(Nbytes);
int i,j,ij;
for(i=0;i<Nx;i++){
for(j=0;j<Ny;j++){
ij = i*Ny+j;
a[ij] = 1.0;
b[ij] = 2.0;
c[ij] = 0.0;
}
}
init(a, b, c);
add(a, b, c);
free(a);
free(b);
free(c);
return 0;
}
void add(float *a, float *b, float *c){
int i,j,ij;
for(i=0;i<Nx;i++){
for(j=0;j<Ny;j++){
576
}
vectoradd2d_1darray.c
GPGPU実践基礎工学
2015/11/25
ベクトル和（2次元並列版）
#define Nx (1024)
#define Ny (1024)
#define Nbytes (Nx*Ny*sizeof(float))
#define NTx (16)
#define NTy (16)
#define NBx (Nx/NTx)
#define NBy (Ny/NTy)
int ij = i*Ny + j;
c[ij] = a[ij] + b[ij];
}
__global__ void init(float *a, float *b,
float *c){
int main(void){
float *a,*b,*c;
dim3 thread(NTx, NTy, 1);
dim3 block(NBx, NBy, 1);
int i = blockIdx.x*blockDim.x + threadIdx.x;
int j = blockIdx.y*blockDim.y + threadIdx.y;
int ij = i*Ny + j;
cudaMalloc( (void **)&a, Nbytes);
cudaMalloc( (void **)&b, Nbytes);
cudaMalloc( (void **)&c, Nbytes);
a[ij] = 1.0;
b[ij] = 2.0;
c[ij] = 0.0;
init<<< block, thread >>>(a,b,c);
add<<< block, thread >>>(a,b,c);
cudaFree(a);
cudaFree(b);
cudaFree(c);
return 0;
}
__global__ void add(float *a, float *b,
float *c){
int i = blockIdx.x*blockDim.x + threadIdx.x;
int j = blockIdx.y*blockDim.y + threadIdx.y;
577
}
GPGPU実践基礎工学
vectoradd2d.cu
2015/11/25
実行結果

2次元版
# CUDA_PROFILE_LOG_VERSION 2.0
# CUDA_DEVICE 0 Tesla M2050
# TIMESTAMPFACTOR fffff5f7a904f250
method,gputime,cputime,occupancy
method=[ _Z4initPfS_S_ ] gputime=[ 735.072 ] cputime=[ 16.000 ] occupancy=[ 1.000 ]
method=[ _Z3addPfS_S_ ] gputime=[ 298.080 ] cputime=[ 7.000 ] occupancy=[ 1.000 ]
カーネル
init
add
578
スレッド数
16×16
16×16
GPGPU実践基礎工学
実行時間[s]
735
298
2015/11/25
実行結果

1次元版(vectoradd.cu)
# CUDA_PROFILE_LOG_VERSION 2.0
# CUDA_DEVICE 0 Tesla M2050
# TIMESTAMPFACTOR fffff5f7dad0bb38
method,gputime,cputime,occupancy
method=[ _Z4initPfS_S_ ] gputime=[ 108.928 ] cputime=[ 15.000 ] occupancy=[ 1.000 ]
method=[ _Z3addPfS_S_ ] gputime=[ 112.992 ] cputime=[ 7.000 ] occupancy=[ 1.000 ]
カーネル
init
add
579
スレッド数
実行時間[s]
109
256(=16×16)
113
256(=16×16)
GPGPU実践基礎工学
2015/11/25
ベクトル和（2次元並列版）の実行結果

実行結果




1次元版と比較して実行時間がかかる
initで7倍
addで3倍
遅くなる要因は？



580
2次元は添字iだけでなくjも計算
添字の計算負荷は軽い
実行時の並列度の指定
GPGPU実践基礎工学
2015/11/25
ベクトル和（2次元並列版）
#define Nx (1024)
#define Ny (1024)
#define Nbytes (Nx*Ny*sizeof(float))
#define NTx (16)
#define NTy (16)
#define NBx (Nx/NTx)
#define NBy (Ny/NTy)
int ij = i+ Nx*j;
c[ij] = a[ij] + b[ij];
}
__global__ void init(float *a, float *b,
float *c){
int main(void){
float *a,*b,*c;
dim3 thread(NTx, NTy, 1);
dim3 block(NBx, NBy, 1);
int i = blockIdx.x*blockDim.x + threadIdx.x;
int j = blockIdx.y*blockDim.y + threadIdx.y;
int ij = i+ Nx*j;
cudaMalloc( (void **)&a, Nbytes);
cudaMalloc( (void **)&b, Nbytes);
cudaMalloc( (void **)&c, Nbytes);
a[ij] = 1.0;
b[ij] = 2.0;
c[ij] = 0.0;
init<<< block, thread >>>(a,b,c);
add<<< block, thread >>>(a,b,c);
cudaFree(a);
cudaFree(b);
cudaFree(c);
return 0;
}
__global__ void add(float *a, float *b,
float *c){
int i = blockIdx.x*blockDim.x + threadIdx.x;
int j = blockIdx.y*blockDim.y + threadIdx.y;
581
}
GPGPU実践基礎工学
vectoradd2d.cu
2015/11/25
実行結果

2次元版
# CUDA_PROFILE_LOG_VERSION 2.0
# CUDA_DEVICE 0 Tesla M2050
# TIMESTAMPFACTOR fffff5f8098b61e0
method,gputime,cputime,occupancy
method=[ _Z4initPfS_S_ ] gputime=[ 230.976 ] cputime=[ 15.000 ] occupancy=[ 1.000 ]
method=[ _Z3addPfS_S_ ] gputime=[ 152.960 ] cputime=[ 6.000 ] occupancy=[ 1.000 ]
カーネル
init
add
582
スレッド数
16×16
16×16
GPGPU実践基礎工学
実行時間[s]
230
153
2015/11/25
2次元的な配列アクセスの優先方向(C言語)
2次元配列の場合
in[][],out[][]
for(i=0;i<Nx;i++)
for(j=0;j<Ny;j++)
out[i][j]=in[i][j];
for(j=0;j<Ny;j++)
for(i=0;i<Nx;i++)
out[i][j]=in[i][j];
i
Ny

j
Nx
583
GPGPU実践基礎工学
2015/11/25
2次元的な配列アクセスの優先方向(C言語)
2次元配列の場合
in[][],out[][]
for(i=0;i<Nx;i++)
for(j=0;j<Ny;j++)
out[i][j]=in[i][j];
for(j=0;j<Ny;j++)
for(i=0;i<Nx;i++)
out[i][j]=in[i][j];
i
Ny

j
Nx
584
GPGPU実践基礎工学
2015/11/25
2次元的な配列アクセスの優先方向(C言語)
2次元配列の1次元配列的表現
in[],out[]
for(i=0;i<Nx;i++)
for(j=0;j<Ny;j++)
out[i][j]=in[i][j];
for(i=0;i<Nx;i++)
for(j=0;j<Ny;j++)
out[i*Ny+j]=
in[i*Ny+j];
i
Ny

j
Nx
585
GPGPU実践基礎工学
2015/11/25
2次元的な配列アクセスの優先方向
CUDAで2次元的に並列化してアクセスする場合
for(i=0;i<Nx;i++)
for(j=0;j<Ny;j++)
out[i*Ny+j]=
in[i*Ny+j];
in[],out[]
i
threadIdx.x
threadIdx.y
i = blockIdx.x*blockDim.x
+ threadIdx.x;
j = blockIdx.y*blockDim.y
+ threadIdx.y; j
out[i+Nx*j]=in[i+Nx*j];
Ny

Nx
586
GPGPU実践基礎工学
2015/11/25
2次元的な配列アクセスの優先方向

C言語



CUDA C


メモリアドレスは2次元目が連続
その関係を維持するように配列を1次元化
メモリアクセスの際はthreadIdx.xが連続なアドレスにアク
セスする方が高速
CUDA Cで2次元的な処理を行うとき

587
1次元配列は1次元目のメモリアドレスが連続となるように扱っ
た方がよい
GPGPU実践基礎工学
2015/11/25
2次元配列の1次元表現（C言語）

C言語は2次元目が連続

その関係を維持するように1次元化
picture[Nx*Ny]
2次元配列の場合
picutre[1][2] = 11

1次元配列の場合
picture[_______] = 11
1*8 + 2
9
10
11
12
Ny=8

1
2
3
4
5
6
7
j
8
メモリが連続な方向
i
Nx=8
picture[i][j]  picture[i*Ny+j]
588
GPGPU実践基礎工学
2015/11/25
2次元配列の1次元表現（CUDA C）

メモリアクセスの際はthreadIdx.xが連続なアドレスに
アクセスする方が高速
picture[Nx*Ny]
i
2次元配列の場合
picutre[1][2] = 11

1次元配列の場合
9
10
11
12
Ny=8

1
2
3
4
5
6
7
j
8
Nx=8
589
GPGPU実践基礎工学
2015/11/25
2次元配列の1次元表現（CUDA C）

threadIdx.xが連続なアドレスにアクセスする方がメモ
リアクセスが高速
block(0,0)
block(1,0)
threadIdx.x
2次元配列の場合
picutre[1][2] = 11

1次元配列の場合
9
10
11
12
threadIdx.y

1
2
3
4
5
6
7
8
メモリが連続な方向
picture[_______] = 11
1 + 8*2
block(0,1)
block(1,1)
picture[i][j]  picture[i+Nx*j]
590
GPGPU実践基礎工学
2015/11/25
想定されるカーネル

1スレッドが2次元配列の1要素を計算

添字 i,j の要素を担当
#define Nx (1920) //水平方向ピクセル数
#define Ny (1080) //垂直方向ピクセル数
__global__ void filter(...){
int i = blockIdx.x*blockDim.x + threadIdx.x;
int j = blockIdx.y*blockDim.y + threadIdx.y;
picture[i + Nx*j] = ...;
}
591
GPGPU実践基礎工学
2015/11/25
画像処理の流れ

変数を宣言し，mallocでメモリを確保



画像pを作成


pをfilteredにコピーしておく
処理を実行


画像を格納する変数
p
処理後の画像を格納する変数 filtered
関数を作成し，pとfilteredを渡す
filteredの内容を画面に出力
592
GPGPU実践基礎工学
2015/11/25
メイン関数（CPU版）
#include<stdio.h>
#include<stdlib.h>
#define WHITE (255)
#define BLACK (0)
#define WIDTH 128
#define HEIGHT WIDTH
#define Nbytes (WIDTH*HEIGHT*sizeof(unsigned char))
void create(unsigned char *);
void copy(unsigned char *,unsigned char *);
void print(unsigned char *);
int main(void){
unsigned char *p= (unsigned char *)malloc(Nbytes);
unsigned char *filtered= (unsigned char *)malloc(Nbytes);
//ここで空間フィルタのカーネルを宣言
create(p);
copy(p,filtered);
//ここで処理を行い，結果をfilteredに格納
593
GPGPU実践基礎工学
imageproc.c
2015/11/25
メイン関数（CPU版）
//画面に各画素の値を表示
print(filtered);
return 0;
}
//画像の内容をコピー
void copy(unsigned char *src, unsigned char *dst){
int i;
for(i=0;i<WIDTH*HEIGHT;i++)
dst[i] = src[i];
}
//画像の内容を画面に出力（gnuplotで表示する用）
void print(unsigned char *p){
int i,j;
for(j=0;j<HEIGHT;j++){
for(i=0;i<WIDTH;i++){
printf("%d %d %d¥n", i, j, p[i+WIDTH*j]);//横位置，縦位置，画像の色情報
}
printf("¥n");//1行分表示したら改行を入れる（gnuplotで必要）
}
imageproc.c
}
594
GPGPU実践基礎工学
2015/11/25
画像の作成
void create(unsigned char *p){
6WIDTH/16
x_origin = 9*WIDTH/16;
y_origin = 9*HEIGHT/16;
for(i=0;i<6*WIDTH/16;i++)
for(j=i;j<6*HEIGHT/16;j++)
p[i+x_origin + WIDTH*(j+y_origin)] = BLACK;
}
HEIGHT
9HEIGHT/16
for(j=0;j<HEIGHT;j++){
for(i=0;i<WIDTH;i++){
p[i+WIDTH*j] = WHITE;
}
}
6HEIGHT/16
int i,j;
int x_origin,y_origin;
9WIDTH/16
WIDTH
595
GPGPU実践基礎工学
2015/11/25
結果の画面表示とファイル保存

print関数が画面に各画素の座標と値を表示
0 0 255
1 0 255
:
126 127 255
127 127 255

画面出力のファイル表示（出力のリダイレクトを利用）


>を利用して保存
$ コマンド > ファイル名
$ ./a.out > pic.txt
$ cat pic.txt
0 0 255
1 0 255
:
126 127 255
127 127 255
596
GPGPU実践基礎工学
2015/11/25
gnuplotによる結果の表示

2次元，3次元データをプロットするアプリケーション



コマンドラインで命令を実行してグラフを描画
関数の描画，ファイルから読み込んだデータの表示が可能
tesla??では正しく動作しないため，grouseで実行する
こと
597
GPGPU実践基礎工学
2015/11/25
gnuplotによる結果の表示

表示可能なファイルフォーマット


データはスペース区切り
3次元表示では列（または行）の区切りとして空白行が必要
0 0 255
1 0 255
...
126 0 255
127 0 255
<‐改行
0 1 255
1 1 255
...
126 1 255
127 1 255
<‐改行
0 2 255
1 2 255
598
GPGPU実践基礎工学
2015/11/25
gnuplotによる結果の表示

実行する命令
#表示設定（画像っぽく表示するための設定）
set xrange[0:127] #x軸の表示範囲を0~127（画像サイズ）に固定
set yrange[0:127]
#y軸の表示範囲を0~127（画像サイズ）に固定
set cbrange[0:255]
#カラーバーの範囲を0~255に固定
set size square
#画像を正方形で表示
set pm3d map
#カラー表示して真上から表示
set palette define(0"black",255"white") #0を黒，255を白
#設定を反映してファイルを読み込み
#（ファイル名はシングルクオートで囲む）
splot 'pic.txt'
plot.txt
599
GPGPU実践基礎工学
2015/11/25
画像処理

比較的簡単な4種類の処理を実装

ネガティブ処理

画像反転


空間フィルタ


水平反転，垂直反転
ぼかし，輪郭抽出
モザイク処理
600
GPGPU実践基礎工学
2015/11/25
ネガティブ処理

色の反転処理


黒(0)を白(255)，白(255)を黒(0)に反転
(255−各画素の値)を新しい画素の値とする
原画像
601
処理画像
GPGPU実践基礎工学
2015/11/25
ネガティブ処理
void negative(unsigned char *p, unsigned char *filtered){
int i,j;
for(j=0;j<HEIGHT;j++){
for(i=0;i<WIDTH;i++){
filtered[i+WIDTH*j] = WHITE ‐ p[i + WIDTH*j];
}
}
}
imageproc.c
602
GPGPU実践基礎工学
2015/11/25
水平反転
画像の幾何的位置を反転


HEIGHT‐1

水平反転 i→(WIDTH‐1)‐i
垂直反転 j→(HEIGHT‐1)‐j 原画像
処理画像
j
0
603
i
WIDTH‐1
GPGPU実践基礎工学
2015/11/25
垂直反転
画像の幾何的位置を反転


HEIGHT‐1

水平反転 i→(WIDTH‐1)‐i
垂直反転 j→(HEIGHT‐1)‐j 原画像
処理画像
j
0
604
i
WIDTH‐1
GPGPU実践基礎工学
2015/11/25
水平反転
void hreflect(unsigned char *p, unsigned char *reflected){
int i,j,iref;
for(j=0;j<HEIGHT;j++){
for(i=0;i<WIDTH;i++){
iref = (WIDTH‐1)‐i; //反転後のx座標値
reflected[iref+WIDTH*j] = p[i + WIDTH*j];
}
}
}
imageproc.c
605
GPGPU実践基礎工学
2015/11/25
垂直反転
void vreflect(unsigned char *p, unsigned char *reflected){
int i,j,jref;
for(j=0;j<HEIGHT;j++){
for(i=0;i<WIDTH;i++){
jref = (HEIGHT‐1)‐j; //反転後のy座標値
reflected[i+WIDTH*jref] = p[i + WIDTH*j];
}
}
}
imageproc.c
606
GPGPU実践基礎工学
2015/11/25
空間フィルタ

ある画素とその周囲の画素を使って処理


処理の仕方を規定したカーネルを定義
カーネルは1次元配列で表現
原画像
607
輪郭抽出
フィルタ
（カーネル）
a
b
c
d
e
f
g
h
i
0
1
0
1
‐4
1
0
1
0
GPGPU実践基礎工学
= b+d‐4e+f+h
2015/11/25
空間フィルタ



カーネルは1次元配列で表現し，引数で渡す
ぼかし（平均フィルタ）
1/9
1/9
1/9
1/9
1/9
1/9
1/9
1/9
1/9
float blur[9] ={1.0f/9.0f,1.0f/9.0f,1.0f/9.0f,
1.0f/9.0f,1.0f/9.0f,1.0f/9.0f,
1.0f/9.0f,1.0f/9.0f,1.0f/9.0f};
輪郭抽出（ラプラシアンフィルタ）
0
608
1
0
1
‐4
1
0
1
0
float laplacian[9] ={ 0.0f, 1.0f, 0.0f,
1.0f,‐4.0f, 1.0f,
0.0f, 1.0f, 0.0f};
GPGPU実践基礎工学
2015/11/25
空間フィルタ
void boxfilter(unsigned char *p, unsigned char *filtered, float *filter){
int i,j
int result = BLACK;
//端では何もしない
for(j=1;j<HEIGHT‐1;j++){
for(i=1;i<WIDTH‐1;i++){
result = filter[0]*p[(i‐1) + WIDTH*(j‐1)]
+filter[1]*p[(i ) + WIDTH*(j‐1)]
+filter[2]*p[(i+1) + WIDTH*(j‐1)]
+filter[3]*p[(i‐1) + WIDTH*(j )]
+filter[4]*p[(i ) + WIDTH*(j )]
+filter[5]*p[(i+1) + WIDTH*(j )]
+filter[6]*p[(i‐1) + WIDTH*(j+1)]
+filter[7]*p[(i ) + WIDTH*(j+1)]
+filter[8]*p[(i+1) + WIDTH*(j+1)];
if(result<BLACK) result = ‐result; //数値が負になっていれば‐1をかける
if(result>WHITE) result = WHITE; //数値が255を超えていれば255に収める
filtered[i+WIDTH*j] = (unsigned char)result;
}
}
imageproc.c
}
609
GPGPU実践基礎工学
2015/11/25
モザイク処理

画像を小さな領域に分け，その領域を全て同じ色にする

領域内の全画素を，領域内の画素の平均値に置き換える
原画像
610
処理画像
GPGPU実践基礎工学
2015/11/25
モザイク
void mosaic(unsigned char *f,unsigned char *filtered, int mosaicSize){
int i,j, isub,jsub;
int average;
for(j=0;j<HEIGHT;j+=mosaicSize){
for(i=0;i<WIDTH;i+=mosaicSize){
//領域内の画素の平均値を計算
average = 0;
for(jsub = 0; jsub<mosaicSize; jsub++){
for(isub = 0; isub<mosaicSize; isub++){
average += f[(i+isub) + WIDTH*(j+jsub)];
}
}
average /= (mosaicSize*mosaicSize);
//領域内の画素を計算した平均値で塗りつぶす
for(jsub = 0; jsub<mosaicSize; jsub++){
for(isub = 0; isub<mosaicSize; isub++){
filtered[(i+isub) + WIDTH*(j+jsub)] = (unsigned char)average;
}
}
}
}
imageproc.c
}
611
GPGPU実践基礎工学
2015/11/25
GPUへの移植の方針

画像の作成はCPUで行い，GPUへ転送

1スレッドが1画素の処理を実行


空間フィルタ


2重のforループを排除し，配列添字とスレッド番号を対応
端の画素の処理ができないので今回は処理しない
モザイク


領域の大きさ（横×縦）は1ブロックと同じとする
領域内の画素の平均値の計算は非常に難易度が高い

612
threadIdx.x=0,threadIdx.y=0のスレッドが平均を計算してメモ
リに書き出す
GPGPU実践基礎工学
2015/11/25
メイン関数（GPU版）
#include<stdio.h>
#include<stdlib.h>
#define WHITE (255)
#define BLACK (0)
#define WIDTH 4096
#define HEIGHT WIDTH
#define Nbytes (WIDTH*HEIGHT*sizeof(unsigned char))
#define THREADX 16
#define THREADY 16
#define BLOCKX (WIDTH/THREADX)
#define BLOCKY (HEIGHT/THREADY)
void create(unsigned char *);
void print(unsigned char *);
int main(void){
unsigned char *p= (unsigned char *)malloc(Nbytes);
unsigned char *dev_p, *dev_filtered;
cudaMalloc( (void **)&dev_p, Nbytes);
cudaMalloc( (void **)&dev_filtered, Nbytes);
613
GPGPU実践基礎工学
imageproc.cu
2015/11/25
メイン関数（GPU版）
//ここで空間フィルタのカーネルを宣言
float *filter;
cudaMalloc( (void **)&filter, sizeof(float)*9);
cudaMemcpy(filter, フィルタのカーネル, sizeof(float)*9, cudaMemcpyHostToDevice);
create(p); //CPUで画像を生成してGPUへ送った後，変数dev_filteredにコピーしておく
cudaMemcpy(dev_p, p, Nbytes, cudaMemcpyHostToDevice);
cudaMemcpy(dev_filtered, dev_p, Nbytes, cudaMemcpyDeviceToDevice);
dim3 block( BLOCKX, BLOCKY, 1), thread(THREADX, THREADY, 1);
//ここで処理を行い，結果をdev_filteredに格納
//dev_filteredの内容をCPUへ送る
unsigned char *filtered = (unsigned char *)malloc(Nbytes);
cudaMemcpy(filtered, dev_filtered, Nbytes, cudaMemcpyDeviceToHost);
free(p);
free(filtered);
cudaFree(dev_p);
cudaFree(dev_filtered);
cudaFree(filter);
return 0;
imageproc.cu
}
614
GPGPU実践基礎工学
2015/11/25
メイン関数（GPU版）
void create(unsigned char *p){
int i,j, x_origin,y_origin;
for(j=0;j<HEIGHT;j++){
for(i=0;i<WIDTH;i++){
p[i+WIDTH*j] = WHITE;
}
}
x_origin = 9*WIDTH/16;
y_origin = 9*HEIGHT/16;
for(i=0;i<6*WIDTH/16;i++){
for(j=i;j<6*HEIGHT/16;j++){
p[i+x_origin + WIDTH*(j+y_origin)] = BLACK;
}
}
}
void print(unsigned char *p){
int i,j;
for(j=0;j<HEIGHT;j++){
for(i=0;i<WIDTH;i++){
printf("%d %d %d¥n", i, j, p[i+WIDTH*j]);
}
printf("¥n");
}
}
615
GPGPU実践基礎工学
imageproc.cu
2015/11/25
ネガティブ処理
__global__ void negative(unsigned char *p, unsigned char *filtered){
int i,j;
i = blockIdx.x*blockDim.x + threadIdx.x; //スレッド数と配列添字の対応
j = blockIdx.y*blockDim.y + threadIdx.y; //
filtered[i+WIDTH*j] = WHITE ‐ p[i + WIDTH*j];
}
imageproc.cu
616
GPGPU実践基礎工学
2015/11/25
水平反転
__global__ void hreflect(unsigned char *p, unsigned char *reflected){
int i,j;
int iref;
i = blockIdx.x*blockDim.x + threadIdx.x;
j = blockIdx.y*blockDim.y + threadIdx.y;
iref = (WIDTH‐1)‐i; //反転後のx座標値
reflected[iref+WIDTH*j] = p[i + WIDTH*j];
}
imageproc.cu
617
GPGPU実践基礎工学
2015/11/25
垂直反転
__global__ void vreflect(unsigned char *p, unsigned char *reflected){
int i,j;
int jref;
i = blockIdx.x*blockDim.x + threadIdx.x;
j = blockIdx.y*blockDim.y + threadIdx.y;
jref = (HEIGHT‐1)‐j; //反転後のy座標値
reflected[i+WIDTH*jref] = p[i + WIDTH*j];
}
imageproc.cu
618
GPGPU実践基礎工学
2015/11/25
空間フィルタ
__global__ void boxfilter(unsigned char *p, unsigned char *filtered, float *filter){
int i,j;
int result = BLACK;
i = blockIdx.x*blockDim.x + threadIdx.x;
j = blockIdx.y*blockDim.y + threadIdx.y;
if(0<i && i<WIDTH‐1 && 0<j && j<HEIGHT‐1){ //端の画素は処理をしないため，ifで処理を分岐
result = filter[0]*p[(i‐1) + WIDTH*(j‐1)]
+filter[1]*p[(i ) + WIDTH*(j‐1)]
+filter[2]*p[(i+1) + WIDTH*(j‐1)]
+filter[3]*p[(i‐1) + WIDTH*(j )]
+filter[4]*p[(i ) + WIDTH*(j )]
+filter[5]*p[(i+1) + WIDTH*(j )]
+filter[6]*p[(i‐1) + WIDTH*(j+1)]
+filter[7]*p[(i ) + WIDTH*(j+1)]
+filter[8]*p[(i+1) + WIDTH*(j+1)];
}
if(result<BLACK) result = ‐result; //数値が負になっていれば‐1をかける
if(result>WHITE) result = WHITE; //数値が255を超えていれば255に収める
filtered[i+WIDTH*j] = (unsigned char)result;
}
imageproc.cu
619
GPGPU実践基礎工学
2015/11/25
モザイク処理
__global__ void mosaic(unsigned char *p, unsigned char *filtered, int mosaicSize){
int i,j, isub,jsub;
int average;
i = blockIdx.x*blockDim.x + threadIdx.x;
j = blockIdx.y*blockDim.y + threadIdx.y;
if(threadIdx.x == 0 && threadIdx.y == 0){//ブロック内の1スレッドのみが処理
//領域内の画素の平均値を計算
average = 0;
for(jsub = 0; jsub<mosaicSize; jsub++){
for(isub = 0; isub<mosaicSize; isub++){
average += p[(i+isub) + WIDTH*(j+jsub)];
}
}
average /= (mosaicSize*mosaicSize);
//領域内の画素を計算した平均値で塗りつぶす
for(jsub = 0; jsub<mosaicSize; jsub++){
for(isub = 0; isub<mosaicSize; isub++){
filtered[(i+isub) + WIDTH*(j+jsub)] = (unsigned char)average;
}
}
}
imageproc.cu
}
620
GPGPU実践基礎工学
2015/11/25
処理時間の比較


画像サイズ
4096×4096
1ブロックあたりのスレッド数 16×16

モザイクのサイズはスレッド数と同じ（16×16）
処理
ネガティブ処理
水平反転
垂直反転
空間フィルタ
モザイク処理
621
処理時間[ms]
CPU
175
187
185
553
260
GPGPU実践基礎工学
GPU
1.17
1.18
1.18
4.13
38.5
補足
ビットマップ画像のフォーマット
ビットマップ画像

画像を画素(pixel)の集合として捉え，各画素の色を3
原色（赤緑青，RGB)で表す

1画素あたり1バイト～4バイト

jpg画像のように画像圧縮を施さないのでデータ量が画
像サイズに比例して増大


624
圧縮形式もあるがあまり利用されない
zip形式等のファイル圧縮を使うと効果的に圧縮可能
GPGPU実践基礎工学
2015/11/25
ビットマップファイルフォーマット

ファイルの情報を表すヘッダ＋画像の各画素情報

ヘッダ


ファイルヘッダ


ファイルの種類，ファイルサイズ等
情報ヘッダ


ファイルヘッダ+情報ヘッダの合計54バイト
画像の幅，高さ，1画素あたりのデータサイズなど
画像の各画素の情報
625
GPGPU実践基礎工学
2015/11/25
ビットマップファイルフォーマット

24bitビットマップファイルのバイナリダンプ
ファイルヘッダ
情報ヘッダ
1画素の情報
画像情報
626
GPGPU実践基礎工学
2015/11/25
ビットマップファイルフォーマット

ファイルヘッダ
内容
サイズ
ファイルの種類
備考
u2
ビットマップの場合はBM
u4
16進数表記（リトルエンディアン，下位
ビットが先に書かれる）
予約領域1
u2
必ず0
予約領域2
u2
必ず0
ファイル先頭から画像までの
オフセット
u4
windowsビットマップでは54
（16進数表記で36）
ファイルサイズ(byte)

エンディアン（endian, 多バイトデータの記述・格納形式）
36 00 03 00
リトルエンディアン 00 03 00 3616=196,662 byte
ビッグエンディアン 36 00 03 0016=905,970,432 byte
627
GPGPU実践基礎工学
2015/11/25
ビットマップファイルフォーマット

情報ヘッダ
内容
備考
byte数
情報ヘッダのサイズ
u4
必ず40（byte, 16進数表記で28）
画像の幅(pixel)
4
16進数表記
画像の高さ(pixel)
4
16進数表記
プレーン数
u2
必ず1
色ビット数（bit/pixel）
u2
1,4,8,24,32（16進数表記）
圧縮形式
u4
0,1,2,3のいずれか
画像データサイズ
u4
16進数表記
水平解像度（dot/m）
4
0の場合もある
垂直解像度（dot/m）
4
0の場合もある
パレット数（使用色数）
u4
0の場合もある
重要な色の数
u4
0の場合もある
628
GPGPU実践基礎工学
2015/11/25
ビットマップファイルフォーマット

画像情報（24bitビットマップ）




3原色の情報をそれぞれ1バイト（unsigned char）で保持
3原色の並び順はBGR
幅は行単位で4バイトの倍数に揃えられる
画素の並びは左下から右方向，下から上方向
BG R
629
GPGPU実践基礎工学
2015/11/25

Download Report