CUDAで連続ウェーブレット変換(2D)

CUDAで2次元連続ウェーブレット変換
ソースコード

CUDAで2次元連続ウェーブレット変換†

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

↑

2次元連続ウェーブレット変換†

2次元連続ウェーブレット変換式は以下．

1次元の時に用いたMexican Hatウェーブレット関数の2次元版は，

↑

CUDA実装†

カーネル呼び出し側のコードは以下．

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のウェーブレット変換値を求めます．

カーネルの実装は以下です．

__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);
	}
}

↑