OpenGL入门学习笔记（二）——简单实现FFT海洋
一、前言
文章不赘述OpenGL的使用入门，使用入门请参考LearnOpenGL CN（https://learnopengl-cn.github.io/）。
文章主要参考：
【1】【学习笔记】Unity 基于GPU FFT海洋的实现-理论篇（https://zhuanlan.zhihu.com/p/95482541）
【2】【学习笔记】Unity 基于GPU FFT海洋的实现-实践篇（https://zhuanlan.zhihu.com/p/96811613）
【3】fft海面模拟（二）（https://zhuanlan.zhihu.com/p/64726720）
【4】海面模拟以及渲染（计算着色器、FFT、Reflection Matrix）（https://blog.csdn.net/xiewenzhao123/article/details/79111004）
【5】一小时学会快速傅里叶变换（Fast Fourier Transform）（https://zhuanlan.zhihu.com/p/31584464）
【6】真实感水体渲染技术总结（https://zhuanlan.zhihu.com/p/95917609?utm_source=qq）
二、FFT海洋
基于快速傅立叶变换（Fast Fourier Transform，FFT）的水体渲染方法是目前广泛使用的一种海洋表面渲染方案。
三、代码
const int N = 512;
const int L = 512;
const float spacing = 10.0;
const int BUTTERFLY_STEPS = std::log2(N);
const int LOCAL_WORK_GROUP_SIZE = 32;
const int GLOBAL_WORK_GROUP_SIZE = N / LOCAL_WORK_GROUP_SIZE;const float G = 9.81;
const float A = 0.001;//phillips谱参数，影响波浪高度
const float PI = 3.141592653589;
const glm::vec2 WindDir = glm::vec2(1.0, 2.0);//风向int drawTest();float GeneratePhillipsSpectrum(glm::vec2 k);
float dispersion(glm::vec2 k);glm::vec2 ComputeGaussianRandom(int x, int y);
glm::vec2 gaussian(int x, int y);
unsigned int wangHash(unsigned int seed);
float rand(unsigned int& rngState);
float DonelanBannerDirectionalSpreading(glm::vec2 k);int drawTest()
{auto window = initWindow();Ogle::Shader* heightShader = new Ogle::Shader("./0_1_GenerateHeightSpectrum.comp");Ogle::Shader* displacementShader = new Ogle::Shader("./0_2_GenerateXZDisplacement.comp");Ogle::Shader* normalShader = new Ogle::Shader("./0_3_GenerationNormalBubbles.comp");Ogle::Shader* fftHorShader = new Ogle::Shader("./0_4_FFTHorizontal.comp");Ogle::Shader* fftVerShader = new Ogle::Shader("./0_5_FFTVertical.comp");Ogle::Shader* drawShader = new Ogle::Shader("./0_DrawTest.vert", "./0_DrawTest.frag");float* vertices = new float[N * N * 5];float sX = -(N - 1) * spacing / 2.0;float sZ = (N - 1) * spacing / 2.0;float vX = 0.0;float vY = 0.0;for (int i = 0; i < N; i++){for (int j = 0; j < N; j++){vertices[i * N * 5 + j * 5] = sX;vertices[i * N * 5 + j * 5 + 1] = 0.0;vertices[i * N * 5 + j * 5 + 2] = sZ;vertices[i * N * 5 + j * 5 + 3] = vX;vertices[i * N * 5 + j * 5 + 4] = vY;sX += spacing;vX += 1.0 / (N - 1);}sX = -(N - 1) * spacing / 2.0;sZ -= spacing;vX = 0.0;vY += 1.0 / (N - 1);}unsigned int* indices = new unsigned int[(N - 1) * (N - 1) * 6];for (int i = 0; i < N - 1; i++){for (int j = 0; j < N - 1; j++){indices[i * (N - 1) * 6 + j * 6] = i * N + j;indices[i * (N - 1) * 6 + j * 6 + 1] = i * N + j + 1;indices[i * (N - 1) * 6 + j * 6 + 2] = (i + 1) * N + j;indices[i * (N - 1) * 6 + j * 6 + 3] = i * N + j + 1;indices[i * (N - 1) * 6 + j * 6 + 4] = (i + 1) * N + j;indices[i * (N - 1) * 6 + j * 6 + 5] = (i + 1) * N + j + 1;}}unsigned int VAO;glGenVertexArrays(1, &VAO);glBindVertexArray(VAO);unsigned int VBO;glGenBuffers(1, &VBO);glBindBuffer(GL_ARRAY_BUFFER, VBO);glBufferData(GL_ARRAY_BUFFER, N * N * 5 * 4, vertices, GL_STATIC_DRAW);unsigned int EBO;glGenBuffers(1, &EBO);glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, EBO);glBufferData(GL_ELEMENT_ARRAY_BUFFER, (N - 1) * (N - 1) * 6 * 4, indices, GL_STATIC_DRAW);glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 5 * sizeof(float), (void*)0);glEnableVertexAttribArray(0);glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, 5 * sizeof(float), (void*)(3 * sizeof(float)));glEnableVertexAttribArray(1);///float *h0Data = new float[N * N * 2];float *h0ConjData = new float[N * N * 2];for (int i = 0; i < N; i++) {for (int j = 0; j < N; j++) {//glm::vec2 k(2.0f * PI * (i - N / 2) / L, 2.0f * PI * (j - N / 2) / L);glm::vec2 k(2.0f * PI * i / N - PI, 2.0f * PI * j / N - PI);glm::vec2 gaussian = ComputeGaussianRandom(i, j);glm::vec2 hTilde0 = gaussian * std::sqrt(std::abs(GeneratePhillipsSpectrum(k) * DonelanBannerDirectionalSpreading(k)) / 2.0f);//gaussian = GenerateGaussianRandom();//有可能需要不同的随机数glm::vec2 hTilde0Conj = gaussian * std::sqrt(std::abs(GeneratePhillipsSpectrum(-k) * DonelanBannerDirectionalSpreading(-k)) / 2.0f);hTilde0Conj.y *= -1.0f;//共轭所以y为负h0Data[i * 2 * N + j * 2] = hTilde0.x;h0Data[i * 2 * N + j * 2 + 1] = hTilde0.y;h0ConjData[i * 2 * N + j * 2] = hTilde0Conj.x;h0ConjData[i * 2 * N + j * 2 + 1] = hTilde0Conj.y;}}//存储H0unsigned int g_textureH0;glGenTextures(1, &g_textureH0);//glActiveTexture(GL_TEXTURE0); glBindTexture(GL_TEXTURE_2D, g_textureH0);glTexImage2D(GL_TEXTURE_2D, 0, GL_RG32F, N, N, 0, GL_RG, GL_FLOAT, h0Data);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);//存储H0Conjunsigned int g_textureH0Conj;glGenTextures(1, &g_textureH0Conj);glBindTexture(GL_TEXTURE_2D, g_textureH0Conj);glTexImage2D(GL_TEXTURE_2D, 0, GL_RG32F, N, N, 0, GL_RG, GL_FLOAT, h0ConjData);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);//存储Htunsigned int g_textureHt;glGenTextures(1, &g_textureHt);glBindTexture(GL_TEXTURE_2D, g_textureHt);glTexImage2D(GL_TEXTURE_2D, 0, GL_RG32F, N, N, 0, GL_RG, GL_FLOAT, 0);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);//存储xz便宜图谱unsigned int g_textureDisplacement[2];glGenTextures(1, &g_textureDisplacement[0]);//xglBindTexture(GL_TEXTURE_2D, g_textureDisplacement[0]);glTexImage2D(GL_TEXTURE_2D, 0, GL_RG32F, N, N, 0, GL_RG, GL_FLOAT, 0);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);glGenTextures(1, &g_textureDisplacement[1]);//zglBindTexture(GL_TEXTURE_2D, g_textureDisplacement[1]);glTexImage2D(GL_TEXTURE_2D, 0, GL_RG32F, N, N, 0, GL_RG, GL_FLOAT, 0);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);//存储IFFT后的xyz偏移结果unsigned int g_textureResult[3];glGenTextures(1, &g_textureResult[0]);//xglBindTexture(GL_TEXTURE_2D, g_textureResult[0]);glTexImage2D(GL_TEXTURE_2D, 0, GL_RG32F, N, N, 0, GL_RG, GL_FLOAT, 0);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);glGenTextures(1, &g_textureResult[1]);//yglBindTexture(GL_TEXTURE_2D, g_textureResult[1]);glTexImage2D(GL_TEXTURE_2D, 0, GL_RG32F, N, N, 0, GL_RG, GL_FLOAT, 0);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);glGenTextures(1, &g_textureResult[2]);//zglBindTexture(GL_TEXTURE_2D, g_textureResult[2]);glTexImage2D(GL_TEXTURE_2D, 0, GL_RG32F, N, N, 0, GL_RG, GL_FLOAT, 0);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);//临时存储unsigned int g_textureTemp;glGenTextures(1, &g_textureTemp);glBindTexture(GL_TEXTURE_2D, g_textureTemp);glTexImage2D(GL_TEXTURE_2D, 0, GL_RG32F, N, N, 0, GL_RG, GL_FLOAT, 0);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);//存储法线向量值unsigned int g_textureNormal;glGenTextures(1, &g_textureNormal);glBindTexture(GL_TEXTURE_2D, g_textureNormal);glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA32F, N, N, 0, GL_RGBA, GL_FLOAT, 0);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);//使用GPU计算FFT需要先把第一次迭代的索引计算好例：N=8时 索引为(0,4) (2,6) (1,5) (3,7)float *butterflyIndices = new float[N];for (int i = 0; i < N / 4; i++){butterflyIndices[i * 2] = i * 2;butterflyIndices[i * 2 + 1] = i * 2 + N / 2;butterflyIndices[i * 2 + N / 2] = i * 2 + 1;butterflyIndices[i * 2 + N / 2 + 1] = i * 2 + N / 2 + 1;}//将计算好的索引值存储到贴图当中unsigned int g_textureIndices;glGenTextures(1, &g_textureIndices);glBindTexture(GL_TEXTURE_1D, g_textureIndices);glTexImage1D(GL_TEXTURE_1D, 0, GL_R32F, N, 0, GL_RED, GL_FLOAT, butterflyIndices);glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);///glBindBuffer(GL_ARRAY_BUFFER, 0);glBindVertexArray(0);glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0);glEnable(GL_DEPTH_TEST);while (!glfwWindowShouldClose(window)){float currentFrame = static_cast<float>(glfwGetTime());deltaTime = currentFrame - lastFrame;lastFrame = currentFrame;float timeValue = glfwGetTime();processInput(window);glClearColor(0.0f, 0.0f, 0.0f, 1.0f);glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);//同时清除深度缓冲区///使用HeightSpectrum计算着色器计算高度glUseProgram(heightShader->id);glBindImageTexture(0, g_textureH0, 0, GL_FALSE, 0, GL_READ_ONLY, GL_RG32F);glBindImageTexture(1, g_textureH0Conj, 0, GL_FALSE, 0, GL_READ_ONLY, GL_RG32F);glBindImageTexture(2, g_textureHt, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RG32F);//传递运行时间glUniform1f(glGetUniformLocation(heightShader->id, "time"), timeValue);//创建一个N*N大小的工作组，即同时计算所有的顶点高度值glDispatchCompute(GLOBAL_WORK_GROUP_SIZE, GLOBAL_WORK_GROUP_SIZE, 1);//确保所有的数据都写入到贴图里了glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);//使用XZDisplacement计算着色器计算xz方向偏移glUseProgram(displacementShader->id);glBindImageTexture(0, g_textureHt, 0, GL_FALSE, 0, GL_READ_ONLY, GL_RG32F);glBindImageTexture(1, g_textureDisplacement[0], 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RG32F);glBindImageTexture(2, g_textureDisplacement[1], 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RG32F);//创建一个N*N大小的工作组，即同时计算所有的顶点高度值glDispatchCompute(GLOBAL_WORK_GROUP_SIZE, GLOBAL_WORK_GROUP_SIZE, 1);//确保所有的数据都写入到贴图里了glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);//使用FFT计算着色器对结果进行逆变换//y//HorizontalglUseProgram(fftHorShader->id);glUniform1i(glGetUniformLocation(fftHorShader->id, "u_steps"), BUTTERFLY_STEPS);//绑定索引贴图、上个计算着色器所得结果的波浪函数贴图、将要存储偏移值的贴图glBindImageTexture(0, g_textureHt, 0, GL_FALSE, 0, GL_READ_ONLY, GL_RG32F);glBindImageTexture(1, g_textureTemp, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RG32F);glBindImageTexture(2, g_textureIndices, 0, GL_FALSE, 0, GL_READ_ONLY, GL_R32F);// 先对每一行进行逆变换glDispatchCompute(1, N, 1);glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);//VerticalglUseProgram(fftVerShader->id);glUniform1i(glGetUniformLocation(fftVerShader->id, "u_steps"), BUTTERFLY_STEPS);//绑定索引贴图、上个计算着色器所得结果的波浪函数贴图、将要存储偏移值的贴图glBindImageTexture(0, g_textureTemp, 0, GL_FALSE, 0, GL_READ_ONLY, GL_RG32F);glBindImageTexture(1, g_textureResult[1], 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RG32F);glBindImageTexture(2, g_textureIndices, 0, GL_FALSE, 0, GL_READ_ONLY, GL_R32F);// 再对每一列进行逆变换glDispatchCompute(N, 1, 1);glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);//使用FFT计算着色器对结果进行逆变换//x//HorizontalglUseProgram(fftHorShader->id);glUniform1i(glGetUniformLocation(fftHorShader->id, "u_steps"), BUTTERFLY_STEPS);//绑定索引贴图、上个计算着色器所得结果的波浪函数贴图、将要存储偏移值的贴图glBindImageTexture(0, g_textureDisplacement[0], 0, GL_FALSE, 0, GL_READ_ONLY, GL_RG32F);glBindImageTexture(1, g_textureTemp, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RG32F);glBindImageTexture(2, g_textureIndices, 0, GL_FALSE, 0, GL_READ_ONLY, GL_R32F);// 先对每一行进行逆变换glDispatchCompute(1, N, 1);glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);//VerticalglUseProgram(fftVerShader->id);glUniform1i(glGetUniformLocation(fftVerShader->id, "u_steps"), BUTTERFLY_STEPS);//绑定索引贴图、上个计算着色器所得结果的波浪函数贴图、将要存储偏移值的贴图glBindImageTexture(0, g_textureTemp, 0, GL_FALSE, 0, GL_READ_ONLY, GL_RG32F);glBindImageTexture(1, g_textureResult[0], 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RG32F);glBindImageTexture(2, g_textureIndices, 0, GL_FALSE, 0, GL_READ_ONLY, GL_R32F);// 再对每一列进行逆变换glDispatchCompute(N, 1, 1);glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);//使用FFT计算着色器对结果进行逆变换//z//HorizontalglUseProgram(fftHorShader->id);glUniform1i(glGetUniformLocation(fftHorShader->id, "u_steps"), BUTTERFLY_STEPS);//绑定索引贴图、上个计算着色器所得结果的波浪函数贴图、将要存储偏移值的贴图glBindImageTexture(0, g_textureDisplacement[1], 0, GL_FALSE, 0, GL_READ_ONLY, GL_RG32F);glBindImageTexture(1, g_textureTemp, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RG32F);glBindImageTexture(2, g_textureIndices, 0, GL_FALSE, 0, GL_READ_ONLY, GL_R32F);// 先对每一行进行逆变换glDispatchCompute(1, N, 1);glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);//VerticalglUseProgram(fftVerShader->id);glUniform1i(glGetUniformLocation(fftVerShader->id, "u_steps"), BUTTERFLY_STEPS);//绑定索引贴图、上个计算着色器所得结果的波浪函数贴图、将要存储偏移值的贴图glBindImageTexture(0, g_textureTemp, 0, GL_FALSE, 0, GL_READ_ONLY, GL_RG32F);glBindImageTexture(1, g_textureResult[2], 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RG32F);glBindImageTexture(2, g_textureIndices, 0, GL_FALSE, 0, GL_READ_ONLY, GL_R32F);// 再对每一列进行逆变换glDispatchCompute(N, 1, 1);glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);//更新法线向量glUseProgram(normalShader->id);glBindImageTexture(0, g_textureResult[0], 0, GL_FALSE, 0, GL_READ_ONLY, GL_RG32F);glBindImageTexture(1, g_textureResult[1], 0, GL_FALSE, 0, GL_READ_ONLY, GL_RG32F);glBindImageTexture(2, g_textureResult[2], 0, GL_FALSE, 0, GL_READ_ONLY, GL_RG32F);glBindImageTexture(3, g_textureNormal, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RGBA32F);glDispatchCompute(GLOBAL_WORK_GROUP_SIZE, GLOBAL_WORK_GROUP_SIZE, 1);glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);/glUseProgram(drawShader->id);glUniform3f(glGetUniformLocation(drawShader->id, "lightColor"), 1.0, 1.0, 1.0);glUniform3f(glGetUniformLocation(drawShader->id, "lightPos"), 100.0, 300.0, 100.0);glUniform3f(glGetUniformLocation(drawShader->id, "viewPos"), camera.Position.x, camera.Position.y, camera.Position.z);glm::mat4 model = glm::mat4(1.0);glm::mat4 view = glm::mat4(1.0);glm::mat4 projection = glm::mat4(1.0);view = camera.GetViewMatrix();projection = glm::perspective(glm::radians(camera.Zoom), 800.0f / 600.0f, 0.1f, 10000.0f);glUniformMatrix4fv(glGetUniformLocation(drawShader->id, "model"), 1, GL_FALSE, glm::value_ptr(model));glUniformMatrix4fv(glGetUniformLocation(drawShader->id, "view"), 1, GL_FALSE, glm::value_ptr(view));glUniformMatrix4fv(glGetUniformLocation(drawShader->id, "projection"), 1, GL_FALSE, glm::value_ptr(projection));///将波浪高度贴图绑定在GL_TEXTURE0glActiveTexture(GL_TEXTURE0);glBindTexture(GL_TEXTURE_2D, g_textureResult[0]);//xglActiveTexture(GL_TEXTURE1);glBindTexture(GL_TEXTURE_2D, g_textureResult[1]);//yglActiveTexture(GL_TEXTURE2);glBindTexture(GL_TEXTURE_2D, g_textureResult[2]);//zglActiveTexture(GL_TEXTURE3);glBindTexture(GL_TEXTURE_2D, g_textureNormal);/glBindVertexArray(VAO);//glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);glDrawElements(GL_TRIANGLES, (N - 1) * (N - 1) * 6, GL_UNSIGNED_INT, (void*)0);glfwSwapBuffers(window);glfwPollEvents();}glDeleteVertexArrays(1, &VAO);glDeleteBuffers(1, &VBO);glDeleteBuffers(1, &EBO);glDeleteProgram(drawShader->id);glfwTerminate();return 0;
}float GeneratePhillipsSpectrum(glm::vec2 k)
{float kLength = glm::length(k);kLength = std::max(0.001f, kLength);//kLength = 1;float kLength2 = kLength * kLength;float kLength4 = kLength2 * kLength2;//float k_dot_w = glm::dot(glm::normalize(k), glm::normalize(WindDir));float k_dot_w2 = k_dot_w * k_dot_w;//float windLength = glm::length(WindDir);float l = windLength * windLength / G;float l2 = l * l;//修正因子（Phillips在|k|较大时收敛性较差，使用最后的exp修正参数来抑制较小波浪）float damping = 0.001;float L2 = l2 * damping * damping;//Phillips谱return A * exp(-1.0f / (kLength2 * l2)) / kLength4 /** k_dot_w2*/ * exp(-kLength2 * L2);//*k_dot_w2为方向拓展，可以替换为Donelan-Banner定向传播*
}float dispersion(glm::vec2 k)
{return std::sqrt(G * glm::length(k));
}///
glm::vec2 ComputeGaussianRandom(int x, int y)
{glm::vec2 g = gaussian(x, y);return g;
}glm::vec2 gaussian(int x, int y)
{//均匀分布随机数unsigned int rngState = wangHash(y * N + x);float x1 = rand(rngState);float x2 = rand(rngState);x1 = std::max(1e-6f, x1);x2 = std::max(1e-6f, x2);//计算两个相互独立的高斯随机数float g1 = sqrt(-2.0f * log(x1)) * cos(2.0f * PI * x2);float g2 = sqrt(-2.0f * log(x1)) * sin(2.0f * PI * x2);return glm::vec2(g1, g2);
}unsigned int wangHash(unsigned int seed)
{auto s = (seed ^ 61) ^ (seed >> 16);s *= 9;s = s ^ (s >> 4);s *= 0x27d4eb2d;s = s ^ (s >> 15);return s;
}float rand(unsigned int& rngState)
{// Xorshift算法rngState ^= (rngState << 13);rngState ^= (rngState >> 17);rngState ^= (rngState << 5);return rngState / 4294967296.0f;;
}//Donelan-Banner方向拓展
float DonelanBannerDirectionalSpreading(glm::vec2 k)
{float betaS;float omegap = 0.855f * G / glm::length(WindDir);float ratio = dispersion(k) / omegap;if (ratio < 0.95f){betaS = 2.61f * std::pow(ratio, 1.3f);}if (ratio >= 0.95f && ratio < 1.6f){betaS = 2.28f * std::pow(ratio, -1.3f);}if (ratio > 1.6f){float epsilon = -0.4f + 0.8393f * std::exp(-0.567f * std::log(ratio * ratio));betaS = std::pow(10, epsilon);}float theta = std::atan2(k.y, k.x) - std::atan2(WindDir.y, WindDir.x);return betaS / std::max(1e-7, 2.0f * std::tanh(betaS * PI) * std::pow(std::cosh(betaS * theta), 2));
}
HeightSpectrum计算着色器：
#version 440 coreconst int LOCAL_WORK_GROUP_SIZE = 32;const int L = 512;
const int N = 512;
const float G = 9.81;
const float PI = 3.141592653589;layout(binding = 0 ,rg32f) uniform image2D u_imageH0;
layout(binding = 1 ,rg32f) uniform image2D u_imageH0Conj;
layout(binding = 2 ,rg32f) uniform image2D u_imageOut;uniform float time;layout(local_size_x = LOCAL_WORK_GROUP_SIZE, local_size_y = LOCAL_WORK_GROUP_SIZE, local_size_z = 1) in;float dispersion(vec2 k);
vec2 ComplexMultiply(vec2 complex_1, vec2 complex_2);void main()
{//通过gl_GlobalInvocationID来得知当前执行单元在全局工作组中的位置ivec2 storePos = ivec2(int(gl_GlobalInvocationID.x), int(gl_GlobalInvocationID.y));//vec2 k = vec2(2.0f * PI * (storePos.x - N/2) / L, 2.0f * PI * (storePos.y - N/2) / L);vec2 k = vec2(2.0f * PI * storePos.x / N - PI, 2.0f * PI * storePos.y / N - PI);//根据位置storePos在贴图中采样得到数据vec2 h0 = imageLoad(u_imageH0, storePos).xy;vec2 h0Conj = imageLoad(u_imageH0Conj, storePos).xy;float omegat = dispersion(k) * time;float c = cos(omegat);float s = sin(omegat);vec2 e1 = vec2(c,s);vec2 e2 = vec2(c,-s);vec2 h1 = ComplexMultiply(h0, e1);vec2 h2 = ComplexMultiply(h0Conj, e2);vec2 HTilde = h1 + h2;//将算出来的高度值存储到贴图当中imageStore(u_imageOut, storePos, vec4(HTilde, 0.0, 0.0));
}float dispersion(vec2 k)
{return sqrt(G * length(k));
}vec2 ComplexMultiply(vec2 complex_1, vec2 complex_2)
{float real = complex_1.x * complex_2.x + complex_1.y * complex_2.y * -1.0;float imag = complex_1.x * complex_2.y + complex_1.y * complex_2.x;return vec2(real, imag);
}
XZDisplacement计算着色器：
#version 440 coreconst int LOCAL_WORK_GROUP_SIZE = 32;const int L = 512;
const int N = 512;
const float PI = 3.141592653589;layout(binding = 0 ,rg32f) uniform image2D u_imageHt;
layout(binding = 1 ,rg32f) uniform image2D u_imageDx;
layout(binding = 2 ,rg32f) uniform image2D u_imageDz;layout(local_size_x = LOCAL_WORK_GROUP_SIZE, local_size_y = LOCAL_WORK_GROUP_SIZE, local_size_z = 1) in;vec2 ComplexMultiply(vec2 complex_1, vec2 complex_2);void main()
{//通过gl_GlobalInvocationID来得知当前执行单元在全局工作组中的位置ivec2 storePos = ivec2(int(gl_GlobalInvocationID.x), int(gl_GlobalInvocationID.y));//vec2 k = vec2(2.0f * PI * (storePos.x - N/2) / L, 2.0f * PI * (storePos.y - N/2) / L);vec2 k = vec2(2.0f * PI * storePos.x / N - PI, 2.0f * PI * storePos.y / N - PI);k /= max(0.001f, length(k));//根据位置storePos在贴图中采样得到数据vec2 ht = imageLoad(u_imageHt, storePos).xy;vec2 xHTilde = ComplexMultiply(vec2(0.0,-k.x), ht);vec2 zHTilde = ComplexMultiply(vec2(0.0,-k.y), ht);//按公式还差e^ik·x不知道为什么没乘//将算出来的高度值存储到贴图当中imageStore(u_imageDx, storePos, vec4(xHTilde, 0.0, 0.0));imageStore(u_imageDz, storePos, vec4(zHTilde, 0.0, 0.0));
}vec2 ComplexMultiply(vec2 complex_1, vec2 complex_2)
{float real = complex_1.x * complex_2.x + complex_1.y * complex_2.y * -1.0;float imag = complex_1.x * complex_2.y + complex_1.y * complex_2.x;return vec2(real, imag);
}
NormalBubbles计算着色器：
#version 440 coreconst int LOCAL_WORK_GROUP_SIZE = 32;const int L = 5110;
const int N = 512;
const float PI = 3.141592653589;layout(binding = 0 ,rg32f) uniform image2D u_imageDx;
layout(binding = 1 ,rg32f) uniform image2D u_imageDy;
layout(binding = 2 ,rg32f) uniform image2D u_imageDz;
layout(binding = 3 ,rgba32f) uniform image2D u_imageNormal;
//layout(binding = 4 ,rgba32f) uniform image2D u_imageBubbles;layout(local_size_x = LOCAL_WORK_GROUP_SIZE, local_size_y = LOCAL_WORK_GROUP_SIZE, local_size_z = 1) in;void main()
{//通过gl_GlobalInvocationID来得知当前执行单元在全局工作组中的位置ivec2 storePos = ivec2(int(gl_GlobalInvocationID.x), int(gl_GlobalInvocationID.y));//计算法线float uintLength = L / (N - 1.0f);//两点间单位长度//获取当前点，周围4个点的uv坐标ivec2 uvX1 = ivec2((storePos.x - 1 + N) % N, storePos.y);ivec2 uvX2 = ivec2((storePos.x + 1 + N) % N, storePos.y);ivec2 uvZ1 = ivec2(storePos.x, (storePos.y - 1 + N) % N);ivec2 uvZ2 = ivec2(storePos.x, (storePos.y + 1 + N) % N);//以当前点为中心，获取周围4个点的偏移值vec3 x1D = vec3(imageLoad(u_imageDx, uvX1).x, imageLoad(u_imageDy, uvX1).x, imageLoad(u_imageDz, uvX1).x);vec3 x2D = vec3(imageLoad(u_imageDx, uvX2).x, imageLoad(u_imageDy, uvX2).x, imageLoad(u_imageDz, uvX2).x);vec3 z1D = vec3(imageLoad(u_imageDx, uvZ1).x, imageLoad(u_imageDy, uvZ1).x, imageLoad(u_imageDz, uvZ1).x);vec3 z2D = vec3(imageLoad(u_imageDx, uvZ2).x, imageLoad(u_imageDy, uvZ2).x, imageLoad(u_imageDz, uvZ2).x);//以当前点为原点，构建周围4个点的坐标vec3 x1 = vec3(x1D.x - uintLength, x1D.yz);vec3 x2 = vec3(x2D.x + uintLength, x2D.yz);vec3 z1 = vec3(z1D.xy, z1D.z - uintLength);vec3 z2 = vec3(z2D.xy, z2D.z + uintLength);//计算两个切向量vec3 tangentX = x2 - x1;vec3 tangentZ = z2 - z1;//计算法线vec3 normal = normalize(cross(tangentZ, tangentX));//计算泡沫vec3 ddx = x2D - x1D;vec3 ddz = z2D - z1D;//雅可比行列式float jacobian = (1.0f + ddx.x) * (1.0f + ddz.z) - ddx.z * ddz.x;//jacobian = saturate(max(0, BubblesThreshold - saturate(jacobian)) * BubblesScale);imageStore(u_imageNormal, storePos, vec4(normal, 0.0));//imageStore(u_imageBubbles, storePos, vec4(jacobian, jacobian, jacobian, 0.0));
}
FFTHorizontal计算着色器：
#version 440 coreconst int LOCAL_WORK_GROUP_SIZE = 32;const int L = 512;
const int N = 512;
const float PI = 3.141592653589;uniform int u_steps;layout (binding = 0, rg32f) uniform image2D u_imageIn;
layout (binding = 1, rg32f) uniform image2D u_imageOut;layout (binding = 2, r32f) uniform image1D u_imageIndices;//如果一个变量被声明为shared，那么它将被保存到特定的位置，从而对同一个本地工作组内的所有计算着色器请求可见,通常访问共享shared变量的性能会远远好于访问图像或者着色器存储缓存（例如主内存）的性能
shared vec2 sharedStore[N];layout (local_size_x = LOCAL_WORK_GROUP_SIZE * 8, local_size_y = 1, local_size_z = 1) in;vec2 ComplexMultiply(vec2 complex_1, vec2 complex_2);
vec2 RootOfUnitVector(int n, int k);void main()
{int xIndex = int(gl_GlobalInvocationID.x);int yIndex = int(gl_GlobalInvocationID.y);int leftStoreIndex = 2 * xIndex;int rightStoreIndex = 2 * xIndex + 1;//读取索引（每一组有两个索引例如（0,4））int leftLoadIndex = int(imageLoad(u_imageIndices, leftStoreIndex).r);int rightLoadIndex = int(imageLoad(u_imageIndices, rightStoreIndex).r);ivec2 leftLoadPos;ivec2 rightLoadPos;leftLoadPos = ivec2(leftLoadIndex, yIndex);rightLoadPos = ivec2(rightLoadIndex, yIndex);ivec2 leftStorePos;ivec2 rightStorePos;leftStorePos = ivec2(leftStoreIndex, yIndex);rightStorePos = ivec2(rightStoreIndex, yIndex);//从贴图中读取数据vec2 leftValue = imageLoad(u_imageIn, leftLoadPos).xy;vec2 rightValue = imageLoad(u_imageIn, rightLoadPos).xy;//放入到共享缓存中sharedStore[leftStoreIndex] = leftValue;sharedStore[rightStoreIndex] = rightValue;//确保所有数据都存储完毕（否则后续逻辑将无法读到所需的数据，即要保证时序）memoryBarrierShared();barrier();int numberButterfliesInSection = 1;int currentSection = xIndex;int currentButterfly = 0;//计算FFTfor (int currentStep = 0; currentStep < u_steps; currentStep++){ //根据位置来获取该组所需的两个索引int leftIndex = currentButterfly + currentSection * numberButterfliesInSection * 2;int rightIndex = currentButterfly + numberButterfliesInSection + currentSection * numberButterfliesInSection * 2;//从共享缓存中获得数据leftValue = sharedStore[leftIndex];rightValue = sharedStore[rightIndex];vec2 currentW = RootOfUnitVector(numberButterfliesInSection * 2, currentButterfly);vec2 multiply;vec2 addition;vec2 subtraction;multiply = ComplexMultiply(currentW, rightValue); addition = leftValue + multiply;subtraction = leftValue - multiply; if(currentStep == u_steps-1){sharedStore[leftIndex] = -addition;sharedStore[rightIndex] = -subtraction;}else{sharedStore[leftIndex] = addition;sharedStore[rightIndex] = subtraction;}//确保所有数据计算并存储完毕 memoryBarrierShared();//根据蝴蝶算法来改变参数 numberButterfliesInSection *= 2;currentSection /= 2;currentButterfly = xIndex % numberButterfliesInSection;//确保所有的计算着色器都计算完毕barrier();}int x1 = leftStoreIndex - N/2;int x2 = rightStoreIndex - N/2;int c1 = abs(x1)%2==0?1:-1;int c2 = abs(x2)%2==0?1:-1;imageStore(u_imageOut, leftStorePos, vec4(c1 * sharedStore[leftStoreIndex], 0.0, 0.0));imageStore(u_imageOut, rightStorePos, vec4(c2 * sharedStore[rightStoreIndex], 0.0, 0.0));
}vec2 ComplexMultiply(vec2 complex_1, vec2 complex_2)
{float real = complex_1.x * complex_2.x + complex_1.y * complex_2.y * -1.0;float imag = complex_1.x * complex_2.y + complex_1.y * complex_2.x;return vec2(real, imag);
}//转换成单位根向量
vec2 RootOfUnitVector(int n, int k)
{vec2 result;result.x = cos(2.0 * PI * float(k) / float(n));result.y = sin(2.0 * PI * float(k) / float(n));return result;
}
FFTVertical计算着色器：
#version 440 coreconst int LOCAL_WORK_GROUP_SIZE = 32;const int L = 512;
const int N = 512;
const float PI = 3.141592653589;uniform int u_steps;layout (binding = 0, rg32f) uniform image2D u_imageIn;
layout (binding = 1, rg32f) uniform image2D u_imageOut;layout (binding = 2, r32f) uniform image1D u_imageIndices;//如果一个变量被声明为shared，那么它将被保存到特定的位置，从而对同一个本地工作组内的所有计算着色器请求可见,通常访问共享shared变量的性能会远远好于访问图像或者着色器存储缓存（例如主内存）的性能
shared vec2 sharedStore[N];layout (local_size_x = 1, local_size_y = LOCAL_WORK_GROUP_SIZE * 8, local_size_z = 1) in;vec2 ComplexMultiply(vec2 complex_1, vec2 complex_2);
vec2 RootOfUnitVector(int n, int k);void main()
{int xIndex = int(gl_GlobalInvocationID.x);int yIndex = int(gl_GlobalInvocationID.y);int leftStoreIndex = 2 * yIndex;int rightStoreIndex = 2 * yIndex + 1;//读取索引（每一组有两个索引例如（0,4））int leftLoadIndex = int(imageLoad(u_imageIndices, leftStoreIndex).r);int rightLoadIndex = int(imageLoad(u_imageIndices, rightStoreIndex).r);ivec2 leftLoadPos;ivec2 rightLoadPos;leftLoadPos = ivec2(xIndex, leftLoadIndex);rightLoadPos = ivec2(xIndex, rightLoadIndex);ivec2 leftStorePos;ivec2 rightStorePos;leftStorePos = ivec2(xIndex, leftStoreIndex);rightStorePos = ivec2(xIndex, rightStoreIndex);//从贴图中读取数据vec2 leftValue = imageLoad(u_imageIn, leftLoadPos).xy;vec2 rightValue = imageLoad(u_imageIn, rightLoadPos).xy;//放入到共享缓存中sharedStore[leftStoreIndex] = leftValue;sharedStore[rightStoreIndex] = rightValue;//确保所有数据都存储完毕（否则后续逻辑将无法读到所需的数据，即要保证时序）memoryBarrierShared();barrier();int numberButterfliesInSection = 1;int currentSection = yIndex;int currentButterfly = 0;//计算FFTfor (int currentStep = 0; currentStep < u_steps; currentStep++){ //根据位置来获取该组所需的两个索引int leftIndex = currentButterfly + currentSection * numberButterfliesInSection * 2;int rightIndex = currentButterfly + numberButterfliesInSection + currentSection * numberButterfliesInSection * 2;//从共享缓存中获得数据leftValue = sharedStore[leftIndex];rightValue = sharedStore[rightIndex];vec2 currentW = RootOfUnitVector(numberButterfliesInSection * 2, currentButterfly);vec2 multiply;vec2 addition;vec2 subtraction;multiply = ComplexMultiply(currentW, rightValue); addition = leftValue + multiply;subtraction = leftValue - multiply; if(currentStep == u_steps-1){sharedStore[leftIndex] = -addition;sharedStore[rightIndex] = -subtraction;}else{sharedStore[leftIndex] = addition;sharedStore[rightIndex] = subtraction;}//确保所有数据计算并存储完毕 memoryBarrierShared();//根据蝴蝶算法来改变参数 numberButterfliesInSection *= 2;currentSection /= 2;currentButterfly = yIndex % numberButterfliesInSection;//确保所有的计算着色器都计算完毕barrier();}int z1 = leftStoreIndex - N/2;int z2 = rightStoreIndex - N/2;int c1 = abs(z1)%2==0?1:-1;int c2 = abs(z2)%2==0?1:-1;imageStore(u_imageOut, leftStorePos, vec4(c1 * sharedStore[leftStoreIndex], 0.0, 0.0));imageStore(u_imageOut, rightStorePos, vec4(c2 * sharedStore[rightStoreIndex], 0.0, 0.0));
}vec2 ComplexMultiply(vec2 complex_1, vec2 complex_2)
{float real = complex_1.x * complex_2.x + complex_1.y * complex_2.y * -1.0;float imag = complex_1.x * complex_2.y + complex_1.y * complex_2.x;return vec2(real, imag);
}//转换成单位根向量
vec2 RootOfUnitVector(int n, int k)
{vec2 result;result.x = cos(2.0 * PI * float(k) / float(n));result.y = sin(2.0 * PI * float(k) / float(n));return result;
}
顶点和片元着色器：
#version 440 corelayout (location = 0) in vec3 aPos;
layout (location = 1) in vec2 aTexCoord;layout(binding = 0) uniform sampler2D u_displacementMapX;
layout(binding = 1) uniform sampler2D u_displacementMapY;
layout(binding = 2) uniform sampler2D u_displacementMapZ;
layout(binding = 3) uniform sampler2D u_normal;out vec2 textureCoord;
out vec3 verNormal;
out vec3 FragPos;uniform mat4 model;
uniform mat4 view;
uniform mat4 projection;void main()
{vec3 displacementX = vec3(texture(u_displacementMapX, aTexCoord).r, 0.0, 0.0);vec3 displacementY = vec3(0.0, texture(u_displacementMapY, aTexCoord).r, 0.0);vec3 displacementZ = vec3(0.0, 0.0, texture(u_displacementMapZ, aTexCoord).r);vec3 normal = texture(u_normal,aTexCoord).rgb;gl_Position = projection*view*model*vec4(aPos+displacementX+displacementY+displacementZ, 1.0);textureCoord = aTexCoord;verNormal = normal;FragPos = (model*vec4(aPos+displacementX+displacementY+displacementZ, 1.0)).xyz;
}
#version 440 corein vec2 textureCoord;
in vec3 verNormal;
in vec3 FragPos;uniform vec3 lightColor;
uniform vec3 lightPos;
uniform vec3 viewPos;out vec4 FragColor;void main()
{//环境光vec3 ambient = lightColor * 0.6;//漫反射vec3 norm = normalize(verNormal);vec3 lightDir = normalize(lightPos - FragPos);float diff = max(dot(norm, lightDir), 0.0);vec3 diffuse = diff * lightColor; //镜面光vec3 viewDir = normalize(viewPos - FragPos);vec3 reflectDir = reflect(-lightDir, norm); float spec = pow(max(dot(viewDir, reflectDir), 0.0), 32);vec3 specular = 0.5 * spec * lightColor;//菲涅尔float R0 = 0.02;float cosTheta = dot(viewDir, norm);float R = R0 + (1 - R0) * pow(1 - cosTheta, 5.0);vec3 fresnel = vec3(R);vec3 re = (ambient + diffuse + specular) * vec4(0.0627, 0.145, 0.25, 1.0).xyz;re += fresnel * lightColor * 0.5;FragColor = vec4(re, 1.0);
}
四、实现效果
实现效果比较基本和简单，波浪效果感觉不算特别理想，并且只添加了基本的光照效果，后续再继续完善。
五、2023年5月14日纠错
1、蝴蝶变换的起始索引计算有误，之前没仔细看，参考8 point蝶形网络单纯当成了将索引奇偶对半分然后相互交错的规律（8 point的情况：0、4、2、6、1、5、3、7），再看了一遍，发现其实是将原索引值bit反序算出的新索引。
图源：参考文章【3】
所以这里新增函数：（随便写写，因为整个索引只算一次就不考虑算法效率了，大佬见谅）
int bitreverse(int bit, int in)
{int* v = new int[bit];int t = in;for (int i = bit - 1; i >= 0; i--){int s = t / 2;v[i] = t % 2;t = s;}t = 0;for (int j = bit - 1; j >= 0; j--){if (v[j] == 1){t += std::pow(2, j);}}delete[] v;return t;
}
以及修改原来的索引计算：
float* butterflyIndices = new float[N];
for (int i = 0; i < N; i++)
{butterflyIndices[i] = bitreverse(BUTTERFLY_STEPS, i);
}
2、修改参数WindDir为(10.0, 20.0)，修改参数A为0.00005，最终效果：
波浪整体都有了较大的起伏，不再像原来那样是一种许多小水波的感觉，接近了参考文章【1】中的效果。
ver.2.0：OpenGL入门学习笔记（一）ver.2.0——简单实现FFT海洋
本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！