2012年09月12日

離散フーリエ変換による海面の生成 その1

ランダムノイズ関係の続きです。

海面の生成には、FFT(高速フーリエ変換)を用いた方法が主流のような気がします。
探し方が悪いのかもしれませんが、日本語の情報があまり見当たらないので調べながら試していきます。
(良いページがありましたら教えてください。)

フーリエ変換での海面生成では、大体のところで参考文献[1]をそのまま使ったり改良していたりするようです。
以下は和訳というわけではなくて適当に書いていますので間違っている可能性があります。

FFTベースの海面の高さh(x, t)は、座標x=(x, z)の位置において、次の式で計算します。

ex001.png

ここで、tが時間、kが2次元ベクトルk=(kx, kz)で、kx=2πn/Lx, kz=2πm/Lzでnとmが-N/2≦n<N/2,-M/2≦m<M/2です。
FFTの処理で、離散点x=(nLx/N, mLz/M)における高さフィールドを生成します。ここで、h~(k,t)がフーリエ変換後の周波数領域でのフーリエ係数で、これをランダムノイズをうまく使って逆フーリエ変換すれば海面のようなものが得られる、ということです。

このh~(k,t)をどうやって求めるかということですが、h~0(k)を求めておいて、あとは時間tでexpの部分を変化させることで求めます。
最初のh~0(k)は、次の式で求めます。

ex002.png

ζrとζiは、ガウス乱数(正規分布に従った乱数)で、Box-Muller法[2]などで生成できます。今回は平均0、分散1の乱数を生成します。
iは虚数単位ですので、計算機で計算する場合は(複素数を扱えない場合は)分けて計算する必要があります。

Ph(k)は、Phillips spectrumというもので、風の向きなどを反映させることができるようです。式は次のようになります。

ex003.png

ここで、L=V^2/gで、Vが風速、gが重力定数(≒9.8)、w^が風の方向ベクトル、Aが定数、|k^・w^|^2がコサイン要素で、この要素によって風の方向に動いているように見えるそうです。
この式ではkとkがあって分かりにくかったのですが、kはkの大きさ(=sqrt(kx*kx+kz*kz))のようです。

上の式では、|k|の大きさが値に大きく影響してしまうので、l<<Lとなる小さな長さを定義して、次の式の値を掛けることで抑えることができます。

ex004.png

h~0(k)が求まれば、後は時間で変化させていけばいいということになります。h~(k,t)は次式で計算できます。

ex005.png

ここで、w(k)は周波数で、以下の式などで計算できます。

ex006.png

以上で求まったh~(k,t)を逆フーリエ変換すれば、波の高さh(x, t)が求まることになります。

正直ここまで読んでも実装の仕方がよくわからなかったので、参考文献[3]のソースコードを参考にさせていただきました。

[3]では、FFTWという高速なFFTのライブラリを使用しているのですが、自分は自分で作ったDFTの関数を使っています。
実装はwater.cでほぼ行われています。

wave関数で海面の生成を行っています。まず最初に、initializeWaves関数で、h~0(k)の計算を行っています。
その中で、Phillips関数でPhillips Spectrumを計算して、FFTを使用するために実部と虚部を持つFFTW用の構造体に、(25)式を少し変更した形で計算してあります。

次に、時間tに依存したh~(k,t)の計算を行います。(26)式ではexpの形式で書かれていますが、オイラーの公式e^iθ=cosθ+isinθを使って、式の前半(plus)と後半(minus)、実部(.re)と虚部(.im)の計算に分けています。
このソースコードでは、FFTを用いる場合には半分のデータでいいので、N/2までのループを回しています。

フーリエ変換後の係数は、実部のみのデータをフーリエ変換した場合、0番目が低周波成分、実部の残りのデータがN/2を中心とした線対称の値、虚部の残りのデータがN/2を中心とした点対称の値をとりますので、FFTWを使用する場合には半分でいいということだと思います。
この辺りのお話は離散フーリエ変換で検索すると出てくるかと思います。

自分ではまだFFTを実装していないので、0〜N-1まで計算しています。
その結果を赤を実数、緑を虚数として画像化したh~0(k)の結果(左)と、ある時刻tのh~(k,t)を画像化した結果(右)が以下のようになりました。

h0k.png h1k.png

これは、恐らく風の方向が原因で角度は違いますが、h~(k,t)に関して[4]の左図と同じような傾向が得られているので、ここまでは合っているかと思います。

で、IDFTを行って高さ情報にしてみると、以下のような拡大すると市松模様が入ったような結果が得られました。

big001.png

FFTで計算しているわけではないのでバグの可能性もあるのですが、参考にしたソースコードでは、下記のようなコメントが書かれており、(-1)^(i+j)された値になっているので戻す必要があると書かれていました。

ex007.png

これを計算すると、下記のように一応それっぽい結果が得られました。

idft001.png

実行結果

時刻tを変化させた場合の結果です。



三角関数の重ね合わせでは滑らかになりすぎるので、Choppy Wavesというもので鋭い波を作るところもあるのですが、その辺は後日にしたいと思います。

[1] Jerry Tessendorf: "Simulating Ocean Water", http://graphics.ucsd.edu/courses/rendering/2005/jdewall/tessendorf.pdf
[2] "Generating Gaussian Random Numbers", http://www.taygeta.com/random/gaussian.html
[3] http://www.people.fas.harvard.edu/~reichard/water/
[4] http://www.edxgraphics.com/2/post/2011/10/simulating-ocean-waves-with-fft-on-gpu.html
posted by シンドラー at 01:38 | Comment(0) | TrackBack(0) | Image Processing | このブログの読者になる | 更新情報をチェックする
この記事へのコメント
コメントを書く
お名前:

メールアドレス:

ホームページアドレス:

コメント:

認証コード: [必須入力]


※画像の中の文字を半角で入力してください。
※ブログオーナーが承認したコメントのみ表示されます。

この記事へのトラックバック