OpenGL入门学习笔记（二）——简单实现FFT海洋_水体渲染opengl

作者：知新_RL | 2024-03-19 23:22:43

踩

水体渲染opengl

一、前言

文章不赘述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;
		}
	}

	//存储H0
	unsigned 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);

	//存储H0Conj
	unsigned 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);

	//存储Ht
	unsigned 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]);//x
	glBindTexture(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]);//z
	glBindTexture(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]);//x
	glBindTexture(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]);//y
	glBindTexture(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]);//z
	glBindTexture(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
		//Horizontal
		glUseProgram(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);

		//Vertical
		glUseProgram(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
		//Horizontal
		glUseProgram(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);

		//Vertical
		glUseProgram(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
		//Horizontal
		glUseProgram(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);

		//Vertical
		glUseProgram(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_TEXTURE0
		glActiveTexture(GL_TEXTURE0);
		glBindTexture(GL_TEXTURE_2D, g_textureResult[0]);//x
		glActiveTexture(GL_TEXTURE1);
		glBindTexture(GL_TEXTURE_2D, g_textureResult[1]);//y
		glActiveTexture(GL_TEXTURE2);
		glBindTexture(GL_TEXTURE_2D, g_textureResult[2]);//z
		glActiveTexture(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));
}
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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507

HeightSpectrum计算着色器：

#version 440 core

const 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);
}
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
51
52
53
54
55
56
57
58
59

XZDisplacement计算着色器：

#version 440 core

const 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);
}
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

NormalBubbles计算着色器：

#version 440 core

const 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));
}
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
51
52
53
54
55
56
57
58
59

FFTHorizontal计算着色器：

#version 440 core

const 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;
 
	//计算FFT
	for (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;
}
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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134

FFTVertical计算着色器：

#version 440 core

const 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;
 
	//计算FFT
	for (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;
}
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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134

顶点和片元着色器：

#version 440 core

layout (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;
}
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

#version 440 core

in 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);
}
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

四、实现效果

实现效果比较基本和简单，波浪效果感觉不算特别理想，并且只添加了基本的光照效果，后续再继续完善。

在这里插入图片描述

五、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;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

以及修改原来的索引计算：

float* butterflyIndices = new float[N];
for (int i = 0; i < N; i++)
{
	butterflyIndices[i] = bitreverse(BUTTERFLY_STEPS, i);
}
1
2
3
4
5

2、修改参数WindDir为(10.0, 20.0)，修改参数A为0.00005，最终效果：
波浪整体都有了较大的起伏，不再像原来那样是一种许多小水波的感觉，接近了参考文章【1】中的效果。

在这里插入图片描述

ver.2.0：OpenGL入门学习笔记（一）ver.2.0——简单实现FFT海洋

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/知新_RL/article/detail/268983