CUDAでウェーブレットノイズ
をテンプレートにして作成
[
トップ
|
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
----
#contents
----
*CUDAでウェーブレットノイズ [#qbde4bba]
コンピュータグラフィクスではノイズとしてPerlinノイズ
((K. Perlin and E. M. Hoffert, "Hypertexture", SIGGRAPH '...
が良く用いられます．
このPerlinノイズの欠点(エイリアスとディテール消失のトレー...
((R. L. Cook and T. DeRose, "Wavelet noise", SIGGRAPH 200...
です．
以下に2次元のPerlinノイズとWaveletノイズとそれぞれのフー...
|&ref(noise_perlin_fft_1.jpg);|&ref(noise_wavelet_fft_1.j...
|>|CENTER:左;Perlinノイズ,右:Waveletノイズ|
ディテール消失の原因となっている高周波領域への漏れがWavel...
WaveletノイズはCPUでも問題ない速度で計算できます．
しかし，乱流などに用いる際には複数のバンドのノイズを組み...
計算量が増大して来ます．そのため，このWaveletノイズをCUDA...
**CUDA実装 [#n5f340ff]
C言語での実装コードはWaveletノイズの論文の最後に掲載され...
Waveletノイズではノイズタイルと呼ばれるものを先に計算し，
実際のノイズの値はこれのBスプライン補間により求めます．
今回はノイズタイルはCPUで先に計算して，デバイスメモリに転...
***ホストの処理 [#fc7f700f]
まず，ホスト側のコードは以下．
#code(C){{
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*)ti...
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_uNoi...
// カーネル実行エラーのチェック
cutilCheckMsg("Kernel execution failed");
cutilSafeCall(cudaThreadSynchronize());
// デバイスからホストへ結果を転送
size = nx*ny*sizeof(float);
cutilSafeCall(cudaMemcpy(hW, dW, size, cudaMemcpyDeviceT...
// デバイスメモリ解放
cutilSafeCall(cudaFree(dW));
}
}}
先にCuSetWaveletTile関数に前計算したノイズタイルを渡して...
この作業はノイズタイルのサイズを変えない限り1度しか実行さ...
CuWNoise2Dには結果を受け取る用の配列と生成するノイズ画像...
ノイズの帯域(ここではlevelが大きいほど細かなノイズになり...
デバイスメモリを渡して，
ノイズ画像サイズと同じ数のスレッドを生成してカーネルを実...
***カーネルとデバイス関数 [#hb698b34]
デバイス側ではBスプラインによるノイズタイルの補間を行いま...
補間のコードは以下．
#code(C){{
__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)基底関数を計算 [...
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-スプライン基底関数を計算するところを，
#code(C){{
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を...
#code(C){{
__global__
void wavelet_noise_2d(float *tile, int tile_n, float *wt,...
{
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など)...
#code(C){{
__global__
void wavelet_multiband_noise_2d(float *tile, int tile_n, ...
{
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上へ投...
float sigma_n = 0.265;
sigma_m = sqrtf(sigma_n*sigma_m);
// 分散が1になるようにノイズを調整
if(sigma_m) result /= sigma_m;
wt[i+j*nx] = result;
}
}
}}
**計算結果 [#kdc8f6e1]
512×512の2Dノイズ画像を生成して，その計算時間を計測しまし...
多帯域はlevel = 3〜8で計算しました．
||単帯域|多帯域|h
|CPU|32|178|
|GPU|9|17|
単位はミリ秒です．複数帯域版の結果画像を以下に示します．
CENTER:&ref(multi_noise_wavelet.jpg);
*ソースコード [#ic9141e3]
Visual Studio 2010用のソースコードを以下に置く．
#ref(wavelet_noise.zip)
-必要ライブラリ : libpng,libjpeg,zlib,CUDA4.0,fftw,freegl...
-Gaussianノイズの生成に http://www.math.sci.hiroshima-u.a...
-比較用にPerlin noiseも実装．Perlin noiseのコードは
http://www.mrl.nyu.edu/~perlin/doc/oscar.html および, htt...
終了行:
----
#contents
----
*CUDAでウェーブレットノイズ [#qbde4bba]
コンピュータグラフィクスではノイズとしてPerlinノイズ
((K. Perlin and E. M. Hoffert, "Hypertexture", SIGGRAPH '...
が良く用いられます．
このPerlinノイズの欠点(エイリアスとディテール消失のトレー...
((R. L. Cook and T. DeRose, "Wavelet noise", SIGGRAPH 200...
です．
以下に2次元のPerlinノイズとWaveletノイズとそれぞれのフー...
|&ref(noise_perlin_fft_1.jpg);|&ref(noise_wavelet_fft_1.j...
|>|CENTER:左;Perlinノイズ,右:Waveletノイズ|
ディテール消失の原因となっている高周波領域への漏れがWavel...
WaveletノイズはCPUでも問題ない速度で計算できます．
しかし，乱流などに用いる際には複数のバンドのノイズを組み...
計算量が増大して来ます．そのため，このWaveletノイズをCUDA...
**CUDA実装 [#n5f340ff]
C言語での実装コードはWaveletノイズの論文の最後に掲載され...
Waveletノイズではノイズタイルと呼ばれるものを先に計算し，
実際のノイズの値はこれのBスプライン補間により求めます．
今回はノイズタイルはCPUで先に計算して，デバイスメモリに転...
***ホストの処理 [#fc7f700f]
まず，ホスト側のコードは以下．
#code(C){{
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*)ti...
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_uNoi...
// カーネル実行エラーのチェック
cutilCheckMsg("Kernel execution failed");
cutilSafeCall(cudaThreadSynchronize());
// デバイスからホストへ結果を転送
size = nx*ny*sizeof(float);
cutilSafeCall(cudaMemcpy(hW, dW, size, cudaMemcpyDeviceT...
// デバイスメモリ解放
cutilSafeCall(cudaFree(dW));
}
}}
先にCuSetWaveletTile関数に前計算したノイズタイルを渡して...
この作業はノイズタイルのサイズを変えない限り1度しか実行さ...
CuWNoise2Dには結果を受け取る用の配列と生成するノイズ画像...
ノイズの帯域(ここではlevelが大きいほど細かなノイズになり...
デバイスメモリを渡して，
ノイズ画像サイズと同じ数のスレッドを生成してカーネルを実...
***カーネルとデバイス関数 [#hb698b34]
デバイス側ではBスプラインによるノイズタイルの補間を行いま...
補間のコードは以下．
#code(C){{
__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)基底関数を計算 [...
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-スプライン基底関数を計算するところを，
#code(C){{
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を...
#code(C){{
__global__
void wavelet_noise_2d(float *tile, int tile_n, float *wt,...
{
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など)...
#code(C){{
__global__
void wavelet_multiband_noise_2d(float *tile, int tile_n, ...
{
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上へ投...
float sigma_n = 0.265;
sigma_m = sqrtf(sigma_n*sigma_m);
// 分散が1になるようにノイズを調整
if(sigma_m) result /= sigma_m;
wt[i+j*nx] = result;
}
}
}}
**計算結果 [#kdc8f6e1]
512×512の2Dノイズ画像を生成して，その計算時間を計測しまし...
多帯域はlevel = 3〜8で計算しました．
||単帯域|多帯域|h
|CPU|32|178|
|GPU|9|17|
単位はミリ秒です．複数帯域版の結果画像を以下に示します．
CENTER:&ref(multi_noise_wavelet.jpg);
*ソースコード [#ic9141e3]
Visual Studio 2010用のソースコードを以下に置く．
#ref(wavelet_noise.zip)
-必要ライブラリ : libpng,libjpeg,zlib,CUDA4.0,fftw,freegl...
-Gaussianノイズの生成に http://www.math.sci.hiroshima-u.a...
-比較用にPerlin noiseも実装．Perlin noiseのコードは
http://www.mrl.nyu.edu/~perlin/doc/oscar.html および, htt...
ページ名:
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
