----
#contents
----

*CUDAで2次元連続ウェーブレット変換 [#b892d87d]
[[CUDAで連続ウェーブレット変換]]では1次元のウェーブレット変換を行いました．
これを2次元に拡張して，画像データのウェーブレット変換ができるようにします．

**2次元連続ウェーブレット変換 [#na06fe42]
2次元連続ウェーブレット変換式は以下．

#ref(eq_2dcwt.jpg,nolink)

1次元の時に用いたMexican Hatウェーブレット関数の2次元版は，
#ref(eq_2dmexicanhat.jpg,nolink)

**CUDA実装 [#wa708b2d]
カーネル呼び出し側のコードは以下．

#code(C){{
float* CuCWT2DTex(float *img0, int nx, int ny, float s)
{
	// 信号データ
	cudaArray *caSrc;
	float *dW, *hW;
	int size;

	// 結果格納用デバイスメモリの確保
	size = nx*ny*sizeof(float);
	cutilSafeCall(cudaMalloc((void**)&dW, size));
	cutilSafeCall(cudaMemset((void*)dW, 0, size));

	// 結果格納用ホストメモリの確保
	hW = new float[nx*ny];
	
	int res = 4;
	if(res != 1){
		res = (uint)ceil(res/s);
	}

	// 画像データを格納するテクスチャの準備
	// CUDA Arrayの確保とホストからのデータ転送
	cudaChannelFormatDesc cdesc = cudaCreateChannelDesc(32, 0, 0, 0, cudaChannelFormatKindFloat);
	size = nx*ny*sizeof(float);
	cutilSafeCall(cudaMallocArray(&caSrc, &cdesc, nx, ny));
	cutilSafeCall(cudaMemcpyToArray(caSrc, 0, 0, img0, size, cudaMemcpyHostToDevice));

	// テクスチャパラメータ
	g_TexImg.addressMode[0] = cudaAddressModeWrap;
	g_TexImg.addressMode[1] = cudaAddressModeWrap;
	g_TexImg.filterMode = cudaFilterModePoint;
	g_TexImg.normalized = true;	// 正規化されたテクスチャ座標でアクセス

	// CUDA Arrayをテクスチャにバインド
	cutilSafeCall(cudaBindTextureToArray(g_TexImg, caSrc, cdesc));

	// カーネル実行
	dim3 block(BLOCK_SIZE, BLOCK_SIZE);
	dim3 grid((nx+block.x-1)/block.x, (ny+block.y-1)/block.y);

	cwt2dtex<<< grid, block >>>(dW, nx, ny, 0, 0, 1, 1, s, res);

	cutilCheckMsg("Kernel execution failed");	// カーネル実行エラーのチェック
	cutilSafeCall(cudaThreadSynchronize());		// 全てのスレッドが終わるのを待つ

	// デバイスからホストへ結果を転送
	size = nx*ny*sizeof(float);
	cutilSafeCall(cudaMemcpy(hW, dW, size, cudaMemcpyDeviceToHost));

	// メモリ解放
	cutilSafeCall(cudaFreeArray(caSrc));
	cutilSafeCall(cudaFree(dW));

	return hW;
}
}}

画像データはテクスチャメモリに載せています．
スレッド数は画像のピクセル数と同じになるようにして，
各スレッドがその座標でのスケールsのウェーブレット変換値を求めます．

カーネルの実装は以下です．
#code(C){{
__device__ 
float MexicanHat2D(float x, float y)
{
	x = x*x;
	y = y*y;
	return MEXICAN_HAT_C*(x+y-2)*exp(-(x+y)/2.0);
}

__global__ 
void cwt2dtex(float *wt, int nx, int ny, float x0, float y0, float dx, float dy, float s, int res)
{
	// MARK:cwt2dtex
	int i = blockIdx.x*blockDim.x+threadIdx.x;
	int j = blockIdx.y*blockDim.y+threadIdx.y;

	if(i < nx && j < ny){
		float x = (x0+i*dx);
		float y = (y0+j*dy);

		float w = 0.0;

		// ウェーブレット関数の有効範囲
		float kx1 = -s*MEXICAN_HAT_R+x;
		float kx2 =  s*MEXICAN_HAT_R+x;
		float ky1 = -s*MEXICAN_HAT_R+y;
		float ky2 =  s*MEXICAN_HAT_R+y;

		kx1 = fmaxf(kx1, 0.0);
		kx2 = fminf(kx2, nx-1.0);
		ky1 = fmaxf(ky1, 0.0);
		ky2 = fminf(ky2, ny-1.0);

		// convolution
		float kstep = 1.0/res;
		float kx, ky;
		for(kx = kx1; kx <= kx2; kx += kstep){
			float Tx = (kx-x)/s;
			for(ky = ky1; ky <= ky2; ky += kstep){
				float Ty = (ky-y)/s;

				w += tex2D(g_TexImg, kx/nx, ky/ny)*MexicanHat2D(Tx,Ty);
			}
		}
		wt[i+j*nx] = w/(sqrt(s)*res*res);
	}
}
}}

**結果 [#ec4df426]
CPU用にも実装し(OpenMPで2スレッド並列で実行)，計算時間を比較します．
CPU側の実装はウェーブレット関数値の前計算も行っています．
入力データとして，512×512のLena画像をグレースケール化して用いました．

|Scale| cuda | cpu |h
|1| 38 | 5242 |
|4| 39 | 5250 |
|8| 139 | 19880 |
|16| 502 | 76298 |
|32| 1767 | 257270 |

単位はミリ秒です．各スケールでの変換結果を以下に示します(256×256に縮小しています)．

|CENTER:|CENTER:|CENTER:|c
|&ref(lena.jpg);|&ref(lena_wt_s1.jpg);|&ref(lena_wt_s4.jpg);|
|&ref(lena_wt_s8.jpg);|&ref(lena_wt_s16.jpg);|&ref(lena_wt_s32.jpg);|

高速化はできたものの，スケール値を上げたときにタイムアウトエラーが出てプログラムが落ちます．
これについては[[CUDAカーネル実行のタイムアウト]]を参照してください．
本来は多数のスレッドで単純な問題を解くのがGPUに適しているので，タイムアウトしないレベルまで問題をさらに分解する方がよいかも知れません．


*ソースコード [#w25653fe]
Visual Studio 2010用のソースコードを以下に置く．

#ref(cu_cwt.zip)

-要libpng,libjpeg,CUDA4.0
-キーボード'1','2','3'でそれぞれCPU,GPU,GPU(テクスチャメモリ使用)を切り替える．その後，'t'キーで再変換．
-'j','k'でスケーリングの変更．
トップ   編集 差分 履歴 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS