快速傅立叶变换(FFT)的海面模拟

在这篇文章中，我们将根据Tessendorf的论文[1]中的方程来实现统计波浪模型，以模拟海洋水。使用快速傅立叶变换，我们将能够实现实时交互的帧速率。以下提供两个截图，首先是简单的光照渲染，之后的是网格模型的渲染。

文章使用了一个高度函数 $h(\overrightarrow{x},t)$ ，其中 $\overrightarrow{x}$ 为坐标(x,z)(我们这里令y为高度)，时间为t(说明我们的高度会随时间变化)。

其中， $\widetilde{h}$ 是我们的频域函数，我们的海面模拟实际则是将频域空间转换为空域空间(即IDFT)。其中：

L为海面网格尺寸，且有：

由于我们要在计算机中利用索引计算，因此我们需要设 $n^{'}$ 为0，1，2，...，N-1， $m^{'}$ 为0，1，2，...，M-1，因此有：

那么我们的n，m便可以写成 $n^{'}$ ， $m^{'}$ 的表达式：

把 $n^{'}$ ， $m^{'}$ 代入 $\overrightarrow{k}$ 得：

我们现在把最初的高度函数改写成 $n^{'}$ ， $m^{'}$ 的函数：

我们现在看一下它的频域函数 $\widetilde{h}$ ，以及关于 $n^{'}$ ， $m^{'}$ 的函数 $\widetilde{h}^{'}$ :

其中k为 $\overrightarrow{k}$ 的模，t为时间(与h函数中是同一个t)。 $\widetilde{h}_{0}^{*}$ 是 $\widetilde{h}_{0}$ 的共轭复数。

g为重力加速度，代入 $n^{'}$ ， $m^{'}$ 的 $w^{'$ 为：

最后就剩下我们的 $\widetilde{h}_{0}$ 函数和其共轭复数 $\widetilde{h}_{0}^{*}$ 了，我们可以把 $\widetilde{h}_{0}$ 也改写为 $n^{'}$ ， $m^{'}$ 的形式：

其中 $\xi _{r}$ ， $\xi _{i}$ 是独立的高斯随机变量， $P_{h}$ 是菲利普频谱有如下形式：

其A是菲利普频谱参数(影响波高)， $L=V^{2}/g$ 为可能出现的最大波浪，V为风速。同理我们可以吧 $P_{h}$ 改写为 $P_{h}^{'}$ ：

现在，我们应该拥有渲染波高场所需的一切，但是我们将评估用于生成断续波的位移矢量和用于在每个顶点进行照明计算的法线矢量。位移矢量由下式给出：

法线 $\overrightarrow{N}$ 则由下列函数得出：

所有的数学运算已经完成了。那我们的FFT如何应用在这里呢，由于我们的IDFT公式为：

DFT公式为：

我们发现它们几乎一样，IDFT就在e指数上少了负号，而且在前面多了系数1/N。根据为另一篇https://blog.csdn.net/qq_39300235/article/details/103453909讲解的FFT方法则可以直接使用，我们的代码：

#ifndef FFT_H
#define FFT_H#include <math.h>
#include "complex.h"class cFFT {private:unsigned int N;unsigned int log_2_N;float pi2;bool invers_FFT;unsigned int *reversed;complex **T;complex *c;protected:public:cFFT(unsigned int N,bool I);//N为傅立叶N值，I为是否为IFFT~cFFT();unsigned int reverse(unsigned int i);complex t(unsigned int x, unsigned int N);void fft(complex* input, int stride, int offset);
};#endif

类实现在：

#include "fft.h"cFFT::cFFT(unsigned int N,bool I) : N(N), invers_FFT(I),reversed(0), T(0), pi2(2 * M_PI) {c = 0;log_2_N = log(N)/log(2);reversed = new unsigned int[N];     // 预反转for (int i = 0; i < N; i++) reversed[i] = reverse(i);int pow2 = 1;T = new complex*[log_2_N];      // 预处理T（所有W值）for (int i = 0; i < log_2_N; i++) {T[i] = new complex[pow2];for (int j = 0; j < pow2; j++)if (invers_FFT)T[i][j] = t(j, pow2 * 2).conj();//复数的模为1时倒数等于共轭elseT[i][j] = t(j, pow2 * 2);pow2 *= 2;}c = new complex[N];
}cFFT::~cFFT() {if (c) delete [] c;if (T) {for (int i = 0; i < log_2_N; i++) if (T[i]) delete [] T[i];delete [] T;}if (reversed) delete [] reversed;
}unsigned int cFFT::reverse(unsigned int i) {unsigned int res = 0;for (int j = 0; j < log_2_N; j++) {//反序if((i>>j)&1) res|=(1<<(log_2_N-j-1));}return res;
}complex cFFT::t(unsigned int x, unsigned int N) {return complex(cos(-1*pi2 * x / N), sin(-1*pi2 * x / N));
}void cFFT::fft(complex* input, int stride, int offset) {//获取反转后的数据放入c数组for (int i = 0; i < N; i++) c[i] = input[reversed[i] * stride + offset];int w_= 0;//执行第一次循环次数for(int l=2;l<=N;l*=2){int m=l/2;for(complex* p=c;p!=c+N;p+=l){for (int i=0; i<m; i++) {//蝴蝶变换（W在构造函数中计算过保存在T中）complex tW=T[w_][i]*p[i+m];p[i+m]=p[i]-tW;p[i]=p[i]+tW;}}w_++;}for (int i = 0; i < N; i++){input[i * stride + offset] = c[i];}
}

FFT和IFFT的差别也就只需在cFFT::t()函数中进行区分即可，而我们的海面模型正好是不需要1/N系数的，因此该系数省略。

我们接下来看我们海面的的代码：

#ifndef FFT_OCEAN_hpp
#define FFT_OCEAN_hpp#include <stdio.h>
#include <GL/glew.h>
#include <GLFW/glfw3.h>
#include <glm/glm.hpp>
#include <glm/gtc/matrix_transform.hpp>
#include <glm/gtc/type_ptr.hpp>
#include "complex.h"
#include "fft.h"
#include "shader.h"
#include <vector>using namespace glm;struct vertex_ocean{GLfloat x,y,z;//vertexGLfloat nx,ny,nz;//normalGLfloat tx,ty; //texcoordsGLfloat a,b,c;//htilde0GLfloat _a,_b,_c;//htilde0mk conjugateGLfloat ox,oy,oz; // original position
};class cOcean {private:                      // flag to render geometry or surfacefloat g;                                // gravity constantint N, Nplus1;                          // dimension -- N should be a power of 2float A;                                // phillips spectrum parameter -- affects heights of wavesvec2 w;                              // wind parameterfloat length;                           // length parametercomplex *h_tilde,                       // for fast fourier transform*h_tilde_slopex, *h_tilde_slopez,*h_tilde_dx, *h_tilde_dz;cFFT *fft;                              // fast fourier transformfloat stride=0.1;vertex_ocean *vertices;                 // vertices for vertex buffer objectunsigned int *indices;                  // indicies for vertex buffer objectunsigned int indices_count;             // number of indices to renderGLuint VAO,vbo_vertices, ebo_indices;       // vertex buffer objectsShader ocean_Shader,nrm_Shader;GLint vertex, normal, texCooed; // attributes and uniformsprotected:public:cOcean(const int N, const float A, const vec2 w, const float length,const std::string& vs,const std::string& fs,const std::string& nvs,const std::string& nfs,const std::string& ngeo);~cOcean();void release();float dispersion(int n_prime, int m_prime);     // deep waterfloat phillips(int n_prime, int m_prime);       // phillips spectrumcomplex hTilde_0(int n_prime, int m_prime);complex hTilde(float t, int n_prime, int m_prime);void evaluateWavesFFT(float t);void render(float t, glm::vec3 view_pos, glm::mat4 Projection, glm::mat4 View, glm::mat4 Model, bool use_fft);};#endif /* FFT_OCEAN_hpp */

#include "FFT_OCEAN.hpp"float uniformRandomVariable() {return (float)rand()/RAND_MAX;
}complex gaussianRandomVariable() {float x1, x2, w;do {x1 = 2.f * uniformRandomVariable() - 1.f;x2 = 2.f * uniformRandomVariable() - 1.f;w = x1 * x1 + x2 * x2;} while ( w >= 1.f );w = sqrt((-2.f * log(w)) / w);return complex(x1 * w, x2 * w);
}cOcean::cOcean(const int N, const float A, const vec2 w, const float length,const std::string& vs,const std::string& fs,const std::string& nvs,const std::string& nfs,const std::string& ngeo) :
g(9.81), N(N), Nplus1(N+1), A(A), w(w), length(length),
vertices(0), indices(0), h_tilde(0), h_tilde_slopex(0), h_tilde_slopez(0), h_tilde_dx(0), h_tilde_dz(0), fft(0), ocean_Shader(vs.c_str(),fs.c_str()),nrm_Shader(nvs.c_str(),nfs.c_str(),ngeo.c_str())
{h_tilde        = new complex[N*N];h_tilde_slopex = new complex[N*N];h_tilde_slopez = new complex[N*N];h_tilde_dx     = new complex[N*N];h_tilde_dz     = new complex[N*N];fft            = new cFFT(N,1);vertices       = new vertex_ocean[Nplus1*Nplus1];indices        = new unsigned int[Nplus1*Nplus1*10];int index;complex htilde0, htilde0mk_conj;for (int m_prime = 0; m_prime < Nplus1; m_prime++) {for (int n_prime = 0; n_prime < Nplus1; n_prime++) {index = m_prime * Nplus1 + n_prime;htilde0        = hTilde_0( n_prime,  m_prime);htilde0mk_conj = hTilde_0(-n_prime, -m_prime).conj();vertices[index].tx  = float(n_prime)/Nplus1;vertices[index].ty  = float(m_prime)/Nplus1;vertices[index].a  = htilde0.a;vertices[index].b  = htilde0.b;vertices[index]._a = htilde0mk_conj.a;vertices[index]._b = htilde0mk_conj.b;vertices[index].ox = vertices[index].x =  (n_prime - N / 2.0f) * length / N;vertices[index].oy = vertices[index].y =  0.0f;vertices[index].oz = vertices[index].z =  (m_prime - N / 2.0f) * length / N;vertices[index].nx = 0.0f;vertices[index].ny = 1.0f;vertices[index].nz = 0.0f;}}indices_count = 0;for (int m_prime = 0; m_prime < N; m_prime++) {for (int n_prime = 0; n_prime < N; n_prime++) {index = m_prime * Nplus1 + n_prime;indices[indices_count++] = index;               // two trianglesindices[indices_count++] = index + Nplus1;indices[indices_count++] = index + Nplus1 + 1;indices[indices_count++] = index;indices[indices_count++] = index + Nplus1 + 1;indices[indices_count++] = index + 1;}}glGenVertexArrays(1,&VAO);glGenBuffers(1, &vbo_vertices);glGenBuffers(1, &ebo_indices);glBindVertexArray(VAO);glBindBuffer(GL_ARRAY_BUFFER, vbo_vertices);glBufferData(GL_ARRAY_BUFFER, sizeof(vertex_ocean)*(Nplus1)*(Nplus1), nullptr, GL_DYNAMIC_DRAW);glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, ebo_indices);glBufferData(GL_ELEMENT_ARRAY_BUFFER, indices_count*sizeof(unsigned int), indices, GL_STATIC_DRAW);glBindVertexArray(0);glBindBuffer(GL_ARRAY_BUFFER,0);glBindBuffer(GL_ELEMENT_ARRAY_BUFFER,0);ocean_Shader.use();ocean_Shader.setFloat("hei_Lim", A);ocean_Shader.close();vertex=ocean_Shader.getAttibLoaction("aPos");normal=ocean_Shader.getAttibLoaction("normal");texCooed=ocean_Shader.getAttibLoaction("texCoords");
}cOcean::~cOcean() {if (h_tilde)        delete [] h_tilde;if (h_tilde_slopex) delete [] h_tilde_slopex;if (h_tilde_slopez) delete [] h_tilde_slopez;if (h_tilde_dx)     delete [] h_tilde_dx;if (h_tilde_dz)     delete [] h_tilde_dz;if (fft)        delete fft;if (vertices)       delete [] vertices;if (indices)        delete [] indices;
}
void cOcean::release() {glDeleteBuffers(1, &ebo_indices);glDeleteBuffers(1, &vbo_vertices);
}float cOcean::dispersion(int n_prime, int m_prime) {float w_0 = 2.0f * M_PI / 200.0f;float kx = M_PI * (2 * n_prime - N) / length;float kz = M_PI * (2 * m_prime - N) / length;return floor(sqrt(g * sqrt(kx * kx + kz * kz)) / w_0) * w_0;
}float cOcean::phillips(int n_prime, int m_prime) {vec2 k(M_PI * (2 * n_prime - N) / length,M_PI * (2 * m_prime - N) / length);float k_length  = glm::length(k);if (k_length < 0.000001) return 0.0;float k_length2 = k_length  * k_length;float k_length4 = k_length2 * k_length2;float k_dot_w   = glm::dot(glm::normalize(k),glm::normalize(w));float k_dot_w2  = k_dot_w * k_dot_w;float w_length  = glm::length(w);float L         = w_length * w_length / g;float L2        = L * L;float damping   = 0.001;float l2        = L2 * damping * damping;return A * exp(-1.0f / (k_length2 * L2)) / k_length4 * k_dot_w2 * exp(-k_length2 * l2);
}complex cOcean::hTilde_0(int n_prime, int m_prime) {complex r = gaussianRandomVariable();return r * sqrt(phillips(n_prime, m_prime) / 2.0f);
}complex cOcean::hTilde(float t, int n_prime, int m_prime) {int index = m_prime * Nplus1 + n_prime;complex htilde0(vertices[index].a,  vertices[index].b);complex htilde0mkconj(vertices[index]._a, vertices[index]._b);float omegat = dispersion(n_prime, m_prime) * t;float cos_ = cos(omegat);float sin_ = sin(omegat);complex c0(cos_,  sin_);complex c1(cos_, -sin_);complex res = htilde0 * c0 + htilde0mkconj * c1;return res;
}void cOcean::render(float t, glm::vec3 view_pos, glm::mat4 Projection, glm::mat4 View, glm::mat4 Model, bool use_fft) {evaluateWavesFFT(t);glBindVertexArray(VAO);glBindBuffer(GL_ARRAY_BUFFER, vbo_vertices);glBufferSubData(GL_ARRAY_BUFFER, 0, sizeof(vertex_ocean) * Nplus1 * Nplus1, vertices);glEnableVertexAttribArray(vertex);glVertexAttribPointer(vertex, 3, GL_FLOAT, GL_FALSE, sizeof(vertex_ocean), 0);glEnableVertexAttribArray(normal);glVertexAttribPointer(normal, 3, GL_FLOAT, GL_FALSE, sizeof(vertex_ocean), (void *)(3*sizeof(GLfloat)));glEnableVertexAttribArray(texCooed);glVertexAttribPointer(texCooed, 2, GL_FLOAT, GL_FALSE, sizeof(vertex_ocean), (void *)(6*sizeof(GLfloat)));ocean_Shader.use();ocean_Shader.setMat4("projection", Projection);ocean_Shader.setMat4("view", View);ocean_Shader.setVec3("viewPos", view_pos);for(int i=0;i<3;i++){for (int j=0;j<3;j++) {Model=glm::translate(glm::mat4(1.0), glm::vec3(i*length*stride,0,j*length*stride));ocean_Shader.setMat4("model", Model);glDrawElements(GL_TRIANGLES, indices_count, GL_UNSIGNED_INT, 0);}}glBindVertexArray(0);glBindBuffer(GL_ARRAY_BUFFER, 0);ocean_Shader.close();
}//FFT只能用于N为2的k次
void cOcean::evaluateWavesFFT(float t) {float kx, kz, len, lambda = -1.0f;int index, index1;for (int m_prime = 0; m_prime < N; m_prime++) {kz = M_PI * (2.0f * m_prime - N) / length;for (int n_prime = 0; n_prime < N; n_prime++) {kx = M_PI*(2 * n_prime - N) / length;len = sqrt(kx * kx + kz * kz);index = m_prime * N + n_prime;h_tilde[index] = hTilde(t, n_prime, m_prime);h_tilde_slopex[index] = h_tilde[index] * complex(0, kx);h_tilde_slopez[index] = h_tilde[index] * complex(0, kz);if (len < 0.000001f) {h_tilde_dx[index]     = complex(0.0f, 0.0f);h_tilde_dz[index]     = complex(0.0f, 0.0f);} else {h_tilde_dx[index]     = h_tilde[index] * complex(0, -kx/len);h_tilde_dz[index]     = h_tilde[index] * complex(0, -kz/len);}}}for (int m_prime = 0; m_prime < N; m_prime++) {fft->fft(h_tilde,  1, m_prime * N);fft->fft(h_tilde_slopex, 1, m_prime * N);fft->fft(h_tilde_slopez,  1, m_prime * N);fft->fft(h_tilde_dx, 1, m_prime * N);fft->fft(h_tilde_dz, 1, m_prime * N);}for (int n_prime = 0; n_prime < N; n_prime++) {fft->fft(h_tilde, N, n_prime);fft->fft(h_tilde_slopex, N, n_prime);fft->fft(h_tilde_slopez, N, n_prime);fft->fft(h_tilde_dx, N, n_prime);fft->fft(h_tilde_dz, N, n_prime);}int sign;float signs[] = { 1.0f, -1.0f };vec3 n;for (int m_prime = 0; m_prime < N; m_prime++) {for (int n_prime = 0; n_prime < N; n_prime++) {index  = m_prime * N + n_prime;        // index into h_tilde..index1 = m_prime * Nplus1 + n_prime;    // index into verticessign = signs[(n_prime + m_prime) & 1];h_tilde[index]     = h_tilde[index] * sign;// heightvertices[index1].y = stride*h_tilde[index].a;// displacementh_tilde_dx[index] = h_tilde_dx[index] * sign;h_tilde_dz[index] = h_tilde_dz[index] * sign;vertices[index1].x = stride*(vertices[index1].ox + h_tilde_dx[index].a * lambda);vertices[index1].z = stride*(vertices[index1].oz + h_tilde_dz[index].a * lambda);// normalh_tilde_slopex[index] = h_tilde_slopex[index] * sign;h_tilde_slopez[index] = h_tilde_slopez[index] * sign;n = glm::normalize(vec3(0.0f - h_tilde_slopex[index].a, 1.0f, 0.0f - h_tilde_slopez[index].a));vertices[index1].nx =  n.x;vertices[index1].ny =  n.y;vertices[index1].nz =  n.z;// for tilingif (n_prime == 0 && m_prime == 0) {vertices[index1 + N + Nplus1 * N].y = stride*h_tilde[index].a;vertices[index1 + N + Nplus1 * N].x = stride*(vertices[index1 + N + Nplus1 * N].ox + h_tilde_dx[index].a * lambda);vertices[index1 + N + Nplus1 * N].z = stride*(vertices[index1 + N + Nplus1 * N].oz + h_tilde_dz[index].a * lambda);vertices[index1 + N + Nplus1 * N].nx =  n.x;vertices[index1 + N + Nplus1 * N].ny =  n.y;vertices[index1 + N + Nplus1 * N].nz =  n.z;}if (n_prime == 0) {vertices[index1 + N].y = stride*h_tilde[index].a;vertices[index1 + N].x = stride*(vertices[index1 + N].ox + h_tilde_dx[index].a * lambda);vertices[index1 + N].z = stride*(vertices[index1 + N].oz + h_tilde_dz[index].a * lambda);vertices[index1 + N].nx =  n.x;vertices[index1 + N].ny =  n.y;vertices[index1 + N].nz =  n.z;}if (m_prime == 0) {vertices[index1 + Nplus1 * N].y = stride*h_tilde[index].a;vertices[index1 + Nplus1 * N].x = stride*(vertices[index1 + Nplus1 * N].ox + h_tilde_dx[index].a * lambda);vertices[index1 + Nplus1 * N].z = stride*(vertices[index1 + Nplus1 * N].oz + h_tilde_dz[index].a * lambda);vertices[index1 + Nplus1 * N].nx =  n.x;vertices[index1 + Nplus1 * N].ny =  n.y;vertices[index1 + Nplus1 * N].nz =  n.z;}}}
}

用OPENGL绘制出的效果则为：

References:
1. Tessendorf, Jerry. Simulating Ocean Water. In SIGGRAPH 2002 Course Notes #9 (Simulating Nature: Realistic and Interactive Techniques), ACM Press.

快速傅立叶变换(FFT)的海面模拟相关推荐

快速傅立叶变换(FFT)算法(原来这就是蝶形变换)
快速傅立叶变换(FFT)算法(原来这就是蝶形变换) 为了实现FFT的海面模拟,不得不先撸个FFT算法实现. 离散傅立叶变换(DFT) 学习FFT之前,首先要先了解什么是DFT,我们都知道傅立叶变换是将 ...
JavaScript实现快速傅立叶变换FFT算法（附完整源码）
JavaScript实现快速傅立叶变换FFT算法(附完整源码) radianToDegree.js完整源代码 ComplexNumber.js完整源代码 bitLength.js完整源代码 fastF ...
如何使用计算机实现fft,快速傅立叶变换(FFT)的计算机实现..doc
快速傅立叶变换(FFT)的计算机实现. 信号与系统课程设计 --FFT的计算机实现快速傅里叶变换(FFT)的计算机实现赖智鹏华中科技大学电气与电子工程学院0809班U200811806 Emai ...
快速傅立叶变换fft_使用快速傅立叶变换fft从气候数据中提取季节性模式
快速傅立叶变换fft Meteorology students hardly experience smooth and expeditious data analysis. When comes t ...
神经网络中快速傅立叶变换(FFT)的梯度传递
最近需要在神经网络中构造复数酉矩阵做权重系数,使用了快速傅立叶变换和反变换. 但是FFT不是theano的现成操作模块(有人写过对应的代码,所以应该会很快加进去了),所以想自己去写梯度传递来彻底搞清楚 ...
第一次邂逅快速傅立叶变换(FFT)
为了毕业设计,我要学习JPEG,还有视频压缩技术,在JPEG的时候,我就被前面的DCT给挡住了,现如今我终于写了一个FFT程序,发了我好长的时间.如果说是因为我的无知,还是什么,我对学习这类有关数学的 ...
【快速傅立叶变换fft数论变换ntt学习小记】
概述 fft(快速傅立叶变换)是用来解决多项式乘法的nlog(n)算法,它的主要思想是先把多项式的多项式表达法转化成若干个二维点对(x,y)(点值),把相同x的y乘起来(计算),最后利用这些点对计算出 ...
Opencv学习笔记 - 使用快速傅立叶变换(FFT)检测图像清晰度
通常的图像清晰度检测大都是计算sobel.拉普拉斯算子的方差,不过大多数时候,拉普拉斯算子方法需要进行大量的手动调整,才能定义图像是否被视为模糊.如果您可以控制照明条件,环境和图像捕获过程,则效果很好 ...
为什么要进行傅立叶变换？傅立叶变换究竟有何意义？如何用Matlab实现快速傅立叶变换？
https://www.douban.com/note/164400821/ 写在最前面:本文是我阅读了多篇相关文章后对它们进行分析重组整合而得,绝大部分内容非我所原创.在此向多位原创作者致敬!!! ...

快速傅立叶变换(FFT)的海面模拟

快速傅立叶变换(FFT)的海面模拟

快速傅立叶变换(FFT)的海面模拟相关推荐

最新文章

热门文章