Simulating Ocean Waves with FFT on GPU
来源：互联网 发布：iphone6s连不上4g网络 编辑：程序博客网 时间：2024/04/29 09:27
10/20/2011
9 Comments
Last year, I wrote a demo that renders realistic ocean surface. That demo consists of three major difficulties, which are: 1. An infinite plane mesh that grows to the horizon 2. The height field for ocean waves that displaces the surface 3. Shading the ocean surface. I looked up a number of literatures and combined approaches from different places to finish thisamazing demo. This article will only discuss the second difficulty, building the height field for waves using Fast Fourier Transformation (FFT).
Jerry Tessendorf's paper <Simulating Ocean Waves> proposed an approach to build height field for ocean waves based on a statistical model. This approach was vastly adopted by the animation and game industry. The ocean in the movie <Titanic> is simulated using this approach. The famous CryEngine also uses it to build the displacement map for water surface.
There are many, many complicated equations in Tessenforf's paper. The first time when I read it, it appeared math-intense. However, most of them are not directly related to building the height field for waves, but are important oceangraphic concept, such as Gerstner Wave. Here, I will try to illustrate it in an as simple as possible way and will only discuss several key equations that's directly related to implementing this algorithm.
This algorithm is based on a statistical model, in which wave height is a random variable of horizontal position and time,h(X,t). It decomposes the wave height field into a set of sinus waves with different amplitudes and phases. The model itself provides us with a method to generate these amplitudes and phases, and we use inverse FFT as a mean to quickly evaluate the sum of these sinus waves.
First, let's take a look at the equation of which we need to perform FFT to evaluate the sum:
There are many, many complicated equations in Tessenforf's paper. The first time when I read it, it appeared math-intense. However, most of them are not directly related to building the height field for waves, but are important oceangraphic concept, such as Gerstner Wave. Here, I will try to illustrate it in an as simple as possible way and will only discuss several key equations that's directly related to implementing this algorithm.
This algorithm is based on a statistical model, in which wave height is a random variable of horizontal position and time,h(X,t). It decomposes the wave height field into a set of sinus waves with different amplitudes and phases. The model itself provides us with a method to generate these amplitudes and phases, and we use inverse FFT as a mean to quickly evaluate the sum of these sinus waves.
First, let's take a look at the equation of which we need to perform FFT to evaluate the sum:
in which h(x, t) represents the height at horizontal position x= (x, z) at time t. Each ~h(k, t) is a complex number that represents a set of amplitude and phase at time t, and the vector k points to the direction of travel of the given wave, with a magnitude k dependent on the length of the wave(λ):
If L is the size of the height field we animate, and N is the resolution of the grid(dimention of FFT), then we can get a discrete range of vector k.
In my demo, I took the Direct Compute Demo made by nvidia as reference and set L = 2000 and N = 512. We have found the discrete domain of k, we now need to know how to generate amplitude and phase for each ~h(k, t). Oceangraphic research has proposed a particular spectrum for wind-driven waves larger than capillary waves in a fully developed sea, the Phillips spectrum, which is defined by the following equation.
This is an empirical spectrum, in which L = (V^2)/g, where V is the wind speed and g is the gravitational constant. w^ is the direction of the wind. A is numeric constant that respresents the amplitude of waves. In my demo, I set V = 600 and A = 1.0. Further more, this spectrum can be further optimized by suppressing waves smaller that a small length l << L, which can be achieved by multiply the original spectrm by
I set l = L / 1000.0f.
We now want to generate a set of time-dependent amplitude and phase with Phillips spectrum. To do that, we first need to create a set of amplitudes and phases at time zero, and then we animate the field. The following equation will finish the first task.
We now want to generate a set of time-dependent amplitude and phase with Phillips spectrum. To do that, we first need to create a set of amplitudes and phases at time zero, and then we animate the field. The following equation will finish the first task.
In the equation above, Zr and Zi are two independent draws from a Gaussian random number generator with mean 0 and standard deviation 1. Note this is important. In the begining there was some bugs in my Gaussian RND, and the result I got was chaos. Now, given time t, we can create a field of the frequency amplitudes:
Where ~h*(k, t)=~h(-k, t). Above is the equation for generating each ~h(k, t) at different k and t. The only thing we don't know about know is the wave frequency w(k). For water waves there is a well-known relationship between these frequencies and the magnitude of the corresponding wavevectors, ki. In deep water, where the bottom may be ignored, that relationship is:
Here g is again the gravitational constant and k is the magnitude of vector k.
Now, we have a set animated amplitudes and phases to which we can perform FFT to evaluate the sum. We can also perform FFT to other equations to producing the choppy vectors as well as slope vectors of the height field, but they are not discussed here.
The computational cost of performing FFT with dimensions of 512x512 is extremely high on CPU. Therefore I implemented it on the pixel shader on the GPU, the parallel computing ability of which vastly enhances the efficiency of this algorithm, and enable it to run in real-time.
FFT is a delicate algorithm. The basic idea can be illustrated by the figure below:
Now, we have a set animated amplitudes and phases to which we can perform FFT to evaluate the sum. We can also perform FFT to other equations to producing the choppy vectors as well as slope vectors of the height field, but they are not discussed here.
The computational cost of performing FFT with dimensions of 512x512 is extremely high on CPU. Therefore I implemented it on the pixel shader on the GPU, the parallel computing ability of which vastly enhances the efficiency of this algorithm, and enable it to run in real-time.
FFT is a delicate algorithm. The basic idea can be illustrated by the figure below:
In the figure above, x is being transformed into X, and all the crosses in the figure means:
Basic Butterfly Operation
Since these crosses look very much alike butterflies, FFT is also called "Butterfly Algorithm".
My implementation is simply to first pre-compute all the indices and weights(W) in the figures above and pack the pre-computed data into a texture. At each stage, these indices and weights are fetched in the pixel shader and used to perform a single butterfly operation. Two Ping-Pong textures are used to write data back and forth between stages. The following is the HLSL code for a basic butterfly operation along the horizontal direction.
My implementation is simply to first pre-compute all the indices and weights(W) in the figures above and pack the pre-computed data into a texture. At each stage, these indices and weights are fetched in the pixel shader and used to perform a single butterfly operation. Two Ping-Pong textures are used to write data back and forth between stages. The following is the HLSL code for a basic butterfly operation along the horizontal direction.
Below are figures of ~h(k, t), the animated amplitudes and phases produced by the Phillips spectrum, and h(x, t), the height field for waves.
~h(k, t) in frequence domain and h(x, t) in spatial domain
The final result and more introduction can be found here.
oceanmesh1.objDownload File
9 Comments
iodiot link
5/22/2012 07:01:56
Simply best explanation. Big thanks ;)
jinjian
3/18/2013 08:15:16
can i get the program?
MegaPixel
10/29/2013 16:04:36
Nice article. I got to the point where I've calculated H(k,t)~. Is it correct that during the real time update I see concentric circles radially leaving the center of the image ? ...
Waterlicious
5/3/2014 18:12:49
Wonderful! This really cleared up a lot of small -yet complicated- things up for me. After reading Tessendorf's paper, Keith Lantz' implementation, and Eric Bruneton's gpu fft code, this gave me the spark I needed!
Daniel
5/4/2014 03:00:17
Has anyone an example how the pre-computed butterfly texture should look like? Would be really nice, cause i think that is the problem in my implementation.
Thomas
3/6/2015 12:39:54
Thx for this explanation. It helped me a lot to understand how this works. But there is one thing, I still don't understand. As far as I know, additions of complex numbers can't be performed if the complex numbers are given in polar coordinates. So what would I have to do to calculate the following equation?
~h(k,t) = ~h0(k) exp(iwt) + ~H0*(-k) exp(-iwt)
Fynn
6/21/2015 04:19:59
just complex multiply ~h0(k) with exp(iwt) and ~h0*(-k) with exp(-iwt) (note: exp(ix) = cos(x) + isin(x)). then complex add the results, thats all.
Ozgur
7/10/2015 03:46:53
Nice article!
But i guess you misuse the term "gravitational constant" which is equal to 6.674e−11 and is notated by Big "G". It should not be confused with "small g" (g), which is the local gravitational field of Earth equal to approx. 9.81
broken link
8/5/2015 19:55:07
I get erro from gl_Position = gl_ProjectionMatrix * fragmentPos
10/20/2011
9 Comments
Last year, I wrote a demo that renders realistic ocean surface. That demo consists of three major difficulties, which are: 1. An infinite plane mesh that grows to the horizon 2. The height field for ocean waves that displaces the surface 3. Shading the ocean surface. I looked up a number of literatures and combined approaches from different places to finish thisamazing demo. This article will only discuss the second difficulty, building the height field for waves using Fast Fourier Transformation (FFT).
Jerry Tessendorf's paper <Simulating Ocean Waves> proposed an approach to build height field for ocean waves based on a statistical model. This approach was vastly adopted by the animation and game industry. The ocean in the movie <Titanic> is simulated using this approach. The famous CryEngine also uses it to build the displacement map for water surface.
There are many, many complicated equations in Tessenforf's paper. The first time when I read it, it appeared math-intense. However, most of them are not directly related to building the height field for waves, but are important oceangraphic concept, such as Gerstner Wave. Here, I will try to illustrate it in an as simple as possible way and will only discuss several key equations that's directly related to implementing this algorithm.
This algorithm is based on a statistical model, in which wave height is a random variable of horizontal position and time,h(X,t). It decomposes the wave height field into a set of sinus waves with different amplitudes and phases. The model itself provides us with a method to generate these amplitudes and phases, and we use inverse FFT as a mean to quickly evaluate the sum of these sinus waves.
First, let's take a look at the equation of which we need to perform FFT to evaluate the sum:
There are many, many complicated equations in Tessenforf's paper. The first time when I read it, it appeared math-intense. However, most of them are not directly related to building the height field for waves, but are important oceangraphic concept, such as Gerstner Wave. Here, I will try to illustrate it in an as simple as possible way and will only discuss several key equations that's directly related to implementing this algorithm.
This algorithm is based on a statistical model, in which wave height is a random variable of horizontal position and time,h(X,t). It decomposes the wave height field into a set of sinus waves with different amplitudes and phases. The model itself provides us with a method to generate these amplitudes and phases, and we use inverse FFT as a mean to quickly evaluate the sum of these sinus waves.
First, let's take a look at the equation of which we need to perform FFT to evaluate the sum:
in which h(x, t) represents the height at horizontal position x= (x, z) at time t. Each ~h(k, t) is a complex number that represents a set of amplitude and phase at time t, and the vector k points to the direction of travel of the given wave, with a magnitude k dependent on the length of the wave(λ):
If L is the size of the height field we animate, and N is the resolution of the grid(dimention of FFT), then we can get a discrete range of vector k.
In my demo, I took the Direct Compute Demo made by nvidia as reference and set L = 2000 and N = 512. We have found the discrete domain of k, we now need to know how to generate amplitude and phase for each ~h(k, t). Oceangraphic research has proposed a particular spectrum for wind-driven waves larger than capillary waves in a fully developed sea, the Phillips spectrum, which is defined by the following equation.
This is an empirical spectrum, in which L = (V^2)/g, where V is the wind speed and g is the gravitational constant. w^ is the direction of the wind. A is numeric constant that respresents the amplitude of waves. In my demo, I set V = 600 and A = 1.0. Further more, this spectrum can be further optimized by suppressing waves smaller that a small length l << L, which can be achieved by multiply the original spectrm by
I set l = L / 1000.0f.
We now want to generate a set of time-dependent amplitude and phase with Phillips spectrum. To do that, we first need to create a set of amplitudes and phases at time zero, and then we animate the field. The following equation will finish the first task.
We now want to generate a set of time-dependent amplitude and phase with Phillips spectrum. To do that, we first need to create a set of amplitudes and phases at time zero, and then we animate the field. The following equation will finish the first task.
In the equation above, Zr and Zi are two independent draws from a Gaussian random number generator with mean 0 and standard deviation 1. Note this is important. In the begining there was some bugs in my Gaussian RND, and the result I got was chaos. Now, given time t, we can create a field of the frequency amplitudes:
Where ~h*(k, t)=~h(-k, t). Above is the equation for generating each ~h(k, t) at different k and t. The only thing we don't know about know is the wave frequency w(k). For water waves there is a well-known relationship between these frequencies and the magnitude of the corresponding wavevectors, ki. In deep water, where the bottom may be ignored, that relationship is:
Here g is again the gravitational constant and k is the magnitude of vector k.
Now, we have a set animated amplitudes and phases to which we can perform FFT to evaluate the sum. We can also perform FFT to other equations to producing the choppy vectors as well as slope vectors of the height field, but they are not discussed here.
The computational cost of performing FFT with dimensions of 512x512 is extremely high on CPU. Therefore I implemented it on the pixel shader on the GPU, the parallel computing ability of which vastly enhances the efficiency of this algorithm, and enable it to run in real-time.
FFT is a delicate algorithm. The basic idea can be illustrated by the figure below:
Now, we have a set animated amplitudes and phases to which we can perform FFT to evaluate the sum. We can also perform FFT to other equations to producing the choppy vectors as well as slope vectors of the height field, but they are not discussed here.
The computational cost of performing FFT with dimensions of 512x512 is extremely high on CPU. Therefore I implemented it on the pixel shader on the GPU, the parallel computing ability of which vastly enhances the efficiency of this algorithm, and enable it to run in real-time.
FFT is a delicate algorithm. The basic idea can be illustrated by the figure below:
In the figure above, x is being transformed into X, and all the crosses in the figure means:
Basic Butterfly Operation
Since these crosses look very much alike butterflies, FFT is also called "Butterfly Algorithm".
My implementation is simply to first pre-compute all the indices and weights(W) in the figures above and pack the pre-computed data into a texture. At each stage, these indices and weights are fetched in the pixel shader and used to perform a single butterfly operation. Two Ping-Pong textures are used to write data back and forth between stages. The following is the HLSL code for a basic butterfly operation along the horizontal direction.
My implementation is simply to first pre-compute all the indices and weights(W) in the figures above and pack the pre-computed data into a texture. At each stage, these indices and weights are fetched in the pixel shader and used to perform a single butterfly operation. Two Ping-Pong textures are used to write data back and forth between stages. The following is the HLSL code for a basic butterfly operation along the horizontal direction.
Below are figures of ~h(k, t), the animated amplitudes and phases produced by the Phillips spectrum, and h(x, t), the height field for waves.
~h(k, t) in frequence domain and h(x, t) in spatial domain
The final result and more introduction can be found here.
oceanmesh1.objDownload File
9 Comments
iodiot link
5/22/2012 07:01:56
Simply best explanation. Big thanks ;)
jinjian
3/18/2013 08:15:16
can i get the program?
MegaPixel
10/29/2013 16:04:36
Nice article. I got to the point where I've calculated H(k,t)~. Is it correct that during the real time update I see concentric circles radially leaving the center of the image ? ...
Waterlicious
5/3/2014 18:12:49
Wonderful! This really cleared up a lot of small -yet complicated- things up for me. After reading Tessendorf's paper, Keith Lantz' implementation, and Eric Bruneton's gpu fft code, this gave me the spark I needed!
Daniel
5/4/2014 03:00:17
Has anyone an example how the pre-computed butterfly texture should look like? Would be really nice, cause i think that is the problem in my implementation.
Thomas
3/6/2015 12:39:54
Thx for this explanation. It helped me a lot to understand how this works. But there is one thing, I still don't understand. As far as I know, additions of complex numbers can't be performed if the complex numbers are given in polar coordinates. So what would I have to do to calculate the following equation?
~h(k,t) = ~h0(k) exp(iwt) + ~H0*(-k) exp(-iwt)
Fynn
6/21/2015 04:19:59
just complex multiply ~h0(k) with exp(iwt) and ~h0*(-k) with exp(-iwt) (note: exp(ix) = cos(x) + isin(x)). then complex add the results, thats all.
Ozgur
7/10/2015 03:46:53
Nice article!
But i guess you misuse the term "gravitational constant" which is equal to 6.674e−11 and is notated by Big "G". It should not be confused with "small g" (g), which is the local gravitational field of Earth equal to approx. 9.81
broken link
8/5/2015 19:55:07
I get erro from gl_Position = gl_ProjectionMatrix * fragmentPos
0 0
- Simulating Ocean Waves with FFT on GPU
- Simulating Ocean Water
- Simulating Ocean Water
- FFT水面加入Choppy Waves
- 读《Simulating Physics with Computers》
- Simulating queueing systems with simmer
- GPU Powered DeepLearning with NVIDIA DIGITS on EC2
- Run rbgirshick/py-faster-rcnn on ubuntu with gpu
- 170620 How to install tensorflow with GPU on linux
- C++ AMP: .Massive Data Parallelism on the GPU with Microsoft's C++ AMP (Accelerated Massive Parallel
- caffe setup on ubuntu 15.04 with CPU and NO GPU and NO Python
- FFT快速卷积GPU加速----pycuda
- AMGCL----AMG ON GPU
- Julia On GPU
- Android 模拟器 GPU ON
- CUDA on NVIDIA GPU
- Using+GPU+with+Cuda
- 编译ffmpeg with gpu
- 潘丽云：魏尔斯特拉斯的复变函数思想分析（2009）（2011-01-14 22:34:30）
- Socket编程实践(8) --套接字IO超时设置方法
- SAR图像处理 MSTAR数据库利用问题
- 计算机网络
- 第一次成功实现的WebSocket代码
- Simulating Ocean Waves with FFT on GPU
- 源泉书签入口
- Activity的生命周期
- SDWebImage异步加载图片及缓存的管理与清理 -- 推荐
- 学习嵌入式过程的一些感想和思考
- Rendering underwater scene with EDXRay
- Makefile常见问题汇总
- SD实践代码
- hdu 1827 Summer Holiday
