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

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

noise_perlin_fft_1.jpgnoise_wavelet_fft_1.jpg
左;Perlinノイズ,右:Waveletノイズ

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

CUDA実装

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

ホストの処理

まず,ホスト側のコードは以下.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
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スプラインによるノイズタイルの補間を行います. 補間のコードは以下.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
__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;
 
        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-スプライン基底関数を計算するところを,

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
        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を呼び出すだけです.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
__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など)を作りたいときは以下のようになります.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
__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;
        }
 
                float sigma_n = 0.265;
        sigma_m = sqrtf(sigma_n*sigma_m);
 
                if(sigma_m) result /= sigma_m;
 
        wt[i+j*nx] = result;
    }
}

計算結果

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

単帯域多帯域
CPU32178
GPU917

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

multi_noise_wavelet.jpg

ソースコード

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


*1 K. Perlin and E. M. Hoffert, "Hypertexture", SIGGRAPH '89.code
*2 R. L. Cook and T. DeRose, "Wavelet noise", SIGGRAPH 2005.

添付ファイル: filenoise_wavelet_fft_1.jpg 1816件 [詳細] filewavelet_noise.zip 1509件 [詳細] filenoise_perlin_fft_1.jpg 1875件 [詳細] filemulti_noise_wavelet.jpg 1747件 [詳細] filenoise_wavelet_fft.jpg 862件 [詳細] filenoise_perlin_fft.jpg 915件 [詳細]

トップ   編集 凍結 差分 履歴 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2022-11-30 (水) 13:48:06