CUDAでウェーブレットノイズ

CUDAでウェーブレットノイズ†

コンピュータグラフィクスではノイズとしてPerlinノイズ *1 が良く用いられます． このPerlinノイズの欠点(エイリアスとディテール消失のトレードオフなど)を解決したのがWaveletノイズ *2 です． 以下に2次元のPerlinノイズとWaveletノイズとそれぞれのフーリエ変換を示します．

左;Perlinノイズ,右:Waveletノイズ

ディテール消失の原因となっている高周波領域への漏れがWaveletノイズではカットされていることが分かります． WaveletノイズはCPUでも問題ない速度で計算できます． しかし，乱流などに用いる際には複数のバンドのノイズを組み合わせたturbulence function等を用いるため， 計算量が増大して来ます．そのため，このWaveletノイズをCUDAで実装してみます．

CUDA実装†

C言語での実装コードはWaveletノイズの論文の最後に掲載されていますのでこれを参考にCUDA上で実装します． Waveletノイズではノイズタイルと呼ばれるものを先に計算し， 実際のノイズの値はこれのBスプライン補間により求めます． 今回はノイズタイルはCPUで先に計算して，デバイスメモリに転送しています．

ホストの処理†

まず，ホスト側のコードは以下．

float *g_dNoiseTile = 0;
uint g_uNoiseTileSize = 0;

//! ウェーブレットノイズタイルのセット
void CuSetWaveletTile(float *tile, int n)
{
	if(!g_dNoiseTile && g_uNoiseTileSize != n){
		// デバイスメモリの確保とホストからの転送
		int size = n*n*sizeof(float);
		cutilSafeCall(cudaMalloc((void**)&g_dNoiseTile, size));
		cutilSafeCall(cudaMemcpy((void*)g_dNoiseTile, (void*)tile, size, cudaMemcpyHostToDevice));

		g_uNoiseTileSize = n;
	}
}

void CuWNoise2D(float *hW, int nx, int ny, int level)
{
	printf("\n[CuWNoise2D]\n");

	if(!g_dNoiseTile) return;

	// 信号データ
	float *dW;
	int size;
	float w = powf(2.0f, (float)level);

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

	if(!hW) hW = new float[nx*ny];

	dim3 block(BLOCK_SIZE, BLOCK_SIZE);
	dim3 grid((nx+block.x-1)/block.x, (ny+block.y-1)/block.y);

	wavelet_noise_2d<<< grid, block >>>(g_dNoiseTile, g_uNoiseTileSize, dW, nx, ny, w);

	// カーネル実行エラーのチェック
	cutilCheckMsg("Kernel execution failed");
	cutilSafeCall(cudaThreadSynchronize());

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

	// デバイスメモリ解放
	cutilSafeCall(cudaFree(dW));
}

先にCuSetWaveletTile関数に前計算したノイズタイルを渡しておきます． この作業はノイズタイルのサイズを変えない限り1度しか実行されません． CuWNoise2Dには結果を受け取る用の配列と生成するノイズ画像のサイズ， ノイズの帯域(ここではlevelが大きいほど細かなノイズになります)を渡します． デバイスメモリを渡して， ノイズ画像サイズと同じ数のスレッドを生成してカーネルを実行します．

カーネルとデバイス関数†

デバイス側ではBスプラインによるノイズタイルの補間を行います． 補間のコードは以下．

__device__
float wnoise2d(float *tile, int tile_n, float p[2])
{
	int f[3], c[3];		// filter, noise coef. indices
	int mid[3];

	float w[2][3], t, result = 0;

	// 2次のB-スプライン(quadratic B-spline)基底関数を計算 [t^2/2, 1/2+t-t^2, (1-t)^2/2]
	for(int k = 0; k < 2; ++k){
		mid[k] = ceil(p[k]-0.5);
		t = mid[k]-(p[k]-0.5);

		w[k][0] = t*t/2; 
		w[k][2] = (1-t)*(1-t)/2; 
		w[k][1] = 1-w[k][0]-w[k][2];
	}

	// ノイズタイルを基底関数の値で重み付け補間する
	for(f[1] = -1; f[1] <= 1; ++f[1]){
		for(f[0] = -1; f[0] <= 1; ++f[0]){
			float weight = 1;
			for(int k = 0; k < 2; ++k){
				c[k] = fmodf(mid[k]+f[k], tile_n);
				weight *= w[k][f[k]+1];
			}
			result += weight*tile[c[1]*tile_n+c[0]];
		}
	}

	return result;
}

タイルとタイルサイズ，座標値を渡します． ちなみにWaveletノイズはBスプラインを用いているのでノイズの導関数の計算がとても楽です． 上記コードのB-スプライン基底関数を計算するところを，

if(k == d){
			w[k][0] = -t; 
			w[k][1] = 2*t-1;
			w[k][2] = 1-t;
		}
		else{
			w[k][0] = t*t/2; 
			w[k][2] = (1-t)*(1-t)/2; 
			w[k][1] = 1-w[k][0]-w[k][2];
		}

とすれば，d = 0でx方向，1でy方向の導関数が得られます．

カーネル関数ではレベルに応じて座標値を求めて，wnoise2dを呼び出すだけです．

__global__
void wavelet_noise_2d(float *tile, int tile_n, float *wt, int nx, int ny, float wpos)
{
	int i = blockIdx.x*blockDim.x+threadIdx.x;
	int j = blockIdx.y*blockDim.y+threadIdx.y;

	if(i < nx && j < ny){
		float p[2];
		p[0] = (float)(i+0.5)*wpos/nx;
		p[1] = (float)(j+0.5)*wpos/ny;

		wt[i+j*nx] = wnoise2d(tile, tile_n, p);
	}
}

上のカーネルではある特定の帯域のノイズを計算しますが， 複数の帯域のノイズを合成したもの(turbulence functionなど)を作りたいときは以下のようになります．

__global__
void wavelet_multiband_noise_2d(float *tile, int tile_n, float *wt, int nx, int ny, int first, int nbands)
{
	int i = blockIdx.x*blockDim.x+threadIdx.x;
	int j = blockIdx.y*blockDim.y+threadIdx.y;

	if(i < nx && j < ny){
		float p[2];
		p[0] = (float)(i+0.5)/nx;
		p[1] = (float)(j+0.5)/ny;

		float result = 0;
		float w = powf(2.0, (float)first);
		float q[2];
		float sigma_m = 0;

		for(int b = 0; b < nbands; ++b){
			q[0] = p[0]*w;
			q[1] = p[1]*w;
			result += wnoise2d(tile, tile_n, q)/w;

			sigma_m += 1.0/w*w;
			w *= 2.0;
		}

		// B-Splineの平均分散σNは2Dで0.265，3Dで0.210，2D上へ投影した3Dノイズで0.296
		float sigma_n = 0.265;
		sigma_m = sqrtf(sigma_n*sigma_m);

		// 分散が1になるようにノイズを調整
		if(sigma_m) result /= sigma_m;

		wt[i+j*nx] = result;
	}
}

計算結果†

512×512の2Dノイズ画像を生成して，その計算時間を計測しました． 多帯域はlevel = 3〜8で計算しました．

単位はミリ秒です．複数帯域版の結果画像を以下に示します．

ソースコード†

Visual Studio 2010用のソースコードを以下に置く．

wavelet_noise.zip

• 必要ライブラリ : libpng,libjpeg,zlib,CUDA4.0,fftw,freeglut,glew,FTGL
• Gaussianノイズの生成に http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html のメルセンヌツイスタのコードを用いさせていただいています．
• 比較用にPerlin noiseも実装．Perlin noiseのコードは http://www.mrl.nyu.edu/~perlin/doc/oscar.html および, http://mrl.nyu.edu/~perlin/noise/ のものを用いさせていただいています．