Traktor: More work on FFT ocean.

code/Terrain/OceanComponent.cpp

@@ -87,6 +87,10 @@ bool OceanComponent::create(resource::IResourceManager* resourceManager, render:
 			return false;
 	}
 
+	stcd.format = render::TfR32G32B32A32F;
+	m_spectrumTexture = renderSystem->createSimpleTexture(stcd, T_FILE_LINE_W);
+	m_evolvedSpectrumTexture = renderSystem->createSimpleTexture(stcd, T_FILE_LINE_W);
+
 	AlignedVector< render::VertexElement > vertexElements;
 	vertexElements.push_back(render::VertexElement(render::DataUsage::Position, render::DtFloat2, offsetof(OceanVertex, pos)));
 	vertexElements.push_back(render::VertexElement(render::DataUsage::Custom, render::DtFloat1, offsetof(OceanVertex, edge)));
@@ -201,6 +205,139 @@ void OceanComponent::setup(
 	const world::WorldRenderView& worldRenderView
 )
 {
+
+	/* Do once begin */
+
+	{
+		Ref< render::RenderPass > rp = new render::RenderPass(L"Ocean compute spectrum");
+		rp->addBuild([=, this](const render::RenderGraph&, render::RenderContext* renderContext) {
+			auto renderBlock = renderContext->allocNamed< render::ComputeRenderBlock >(L"Ocean spectrum");
+
+			const render::Shader::Permutation perm(render::getParameterHandle(L"Spectrum"));
+
+			renderBlock->program = m_shaderWave->getProgram(perm).program;
+			renderBlock->workSize[0] = 512;
+			renderBlock->workSize[1] = 512;
+			renderBlock->workSize[2] = 1;
+
+			renderBlock->programParams = renderContext->alloc< render::ProgramParameters >();
+			renderBlock->programParams->beginParameters(renderContext);
+			renderBlock->programParams->setFloatParameter(s_handleWorld_Time, worldRenderView.getTime());
+			renderBlock->programParams->setFloatParameter(s_handleOcean_TileIndex, 0);
+			renderBlock->programParams->setImageViewParameter(s_handleOcean_WaveTexture, m_spectrumTexture, 0);
+			renderBlock->programParams->endParameters(renderContext);
+
+			renderContext->compute(renderBlock);
+			renderContext->compute< render::BarrierRenderBlock >(render::Stage::Compute, render::Stage::Vertex, nullptr, 0);
+		});
+		context.getRenderGraph().addPass(rp);
+	}
+
+	{
+		Ref< render::RenderPass > rp = new render::RenderPass(L"Ocean compute spectrum pack conjugate");
+		rp->addBuild([=, this](const render::RenderGraph&, render::RenderContext* renderContext) {
+			auto renderBlock = renderContext->allocNamed< render::ComputeRenderBlock >(L"Ocean spectrum pack conjugate");
+
+			const render::Shader::Permutation perm(render::getParameterHandle(L"SpectrumPackConjugate"));
+
+			renderBlock->program = m_shaderWave->getProgram(perm).program;
+			renderBlock->workSize[0] = 512;
+			renderBlock->workSize[1] = 512;
+			renderBlock->workSize[2] = 1;
+
+			renderBlock->programParams = renderContext->alloc< render::ProgramParameters >();
+			renderBlock->programParams->beginParameters(renderContext);
+			renderBlock->programParams->setFloatParameter(s_handleWorld_Time, worldRenderView.getTime());
+			renderBlock->programParams->setFloatParameter(s_handleOcean_TileIndex, 0);
+			renderBlock->programParams->setImageViewParameter(s_handleOcean_WaveTexture, m_spectrumTexture, 0);
+			renderBlock->programParams->endParameters(renderContext);
+
+			renderContext->compute(renderBlock);
+			renderContext->compute< render::BarrierRenderBlock >(render::Stage::Compute, render::Stage::Vertex, nullptr, 0);
+		});
+		context.getRenderGraph().addPass(rp);
+	}
+
+	/* Do once end */
+
+	{
+		Ref< render::RenderPass > rp = new render::RenderPass(L"Ocean compute spectrum evolve");
+		rp->addBuild([=, this](const render::RenderGraph&, render::RenderContext* renderContext) {
+			auto renderBlock = renderContext->allocNamed< render::ComputeRenderBlock >(L"Ocean spectrum evolve");
+
+			const render::Shader::Permutation perm(render::getParameterHandle(L"SpectrumEvolve"));
+
+			renderBlock->program = m_shaderWave->getProgram(perm).program;
+			renderBlock->workSize[0] = 512;
+			renderBlock->workSize[1] = 512;
+			renderBlock->workSize[2] = 1;
+
+			renderBlock->programParams = renderContext->alloc< render::ProgramParameters >();
+			renderBlock->programParams->beginParameters(renderContext);
+			renderBlock->programParams->setFloatParameter(s_handleWorld_Time, worldRenderView.getTime());
+			renderBlock->programParams->setFloatParameter(s_handleOcean_TileIndex, 0);
+			renderBlock->programParams->setImageViewParameter(s_handleOcean_WaveTexture, m_spectrumTexture, 0);
+			renderBlock->programParams->setImageViewParameter(s_handleOcean_WaveTexture0, m_evolvedSpectrumTexture, 0);
+			renderBlock->programParams->endParameters(renderContext);
+
+			renderContext->compute(renderBlock);
+			renderContext->compute< render::BarrierRenderBlock >(render::Stage::Compute, render::Stage::Vertex, nullptr, 0);
+		});
+		context.getRenderGraph().addPass(rp);
+	}
+
+	{
+		Ref< render::RenderPass > rp = new render::RenderPass(L"Ocean compute inverse FFT X");
+		rp->addBuild([=, this](const render::RenderGraph&, render::RenderContext* renderContext) {
+			auto renderBlock = renderContext->allocNamed< render::ComputeRenderBlock >(L"Ocean inverse FFT X");
+
+			const render::Shader::Permutation perm(render::getParameterHandle(L"InverseFFT_X"));
+
+			renderBlock->program = m_shaderWave->getProgram(perm).program;
+			renderBlock->workSize[0] = 512;
+			renderBlock->workSize[1] = 512;
+			renderBlock->workSize[2] = 1;
+
+			renderBlock->programParams = renderContext->alloc< render::ProgramParameters >();
+			renderBlock->programParams->beginParameters(renderContext);
+			renderBlock->programParams->setFloatParameter(s_handleWorld_Time, worldRenderView.getTime());
+			renderBlock->programParams->setFloatParameter(s_handleOcean_TileIndex, 0);
+			renderBlock->programParams->setImageViewParameter(s_handleOcean_WaveTexture, m_evolvedSpectrumTexture, 0);
+			renderBlock->programParams->setImageViewParameter(s_handleOcean_WaveTexture0, m_evolvedSpectrumTexture, 0);
+			renderBlock->programParams->endParameters(renderContext);
+
+			renderContext->compute(renderBlock);
+			renderContext->compute< render::BarrierRenderBlock >(render::Stage::Compute, render::Stage::Vertex, nullptr, 0);
+		});
+		context.getRenderGraph().addPass(rp);
+	}
+
+	{
+		Ref< render::RenderPass > rp = new render::RenderPass(L"Ocean compute inverse FFT Y");
+		rp->addBuild([=, this](const render::RenderGraph&, render::RenderContext* renderContext) {
+			auto renderBlock = renderContext->allocNamed< render::ComputeRenderBlock >(L"Ocean inverse FFT Y");
+
+			const render::Shader::Permutation perm(render::getParameterHandle(L"InverseFFT_Y"));
+
+			renderBlock->program = m_shaderWave->getProgram(perm).program;
+			renderBlock->workSize[0] = 512;
+			renderBlock->workSize[1] = 512;
+			renderBlock->workSize[2] = 1;
+
+			renderBlock->programParams = renderContext->alloc< render::ProgramParameters >();
+			renderBlock->programParams->beginParameters(renderContext);
+			renderBlock->programParams->setFloatParameter(s_handleWorld_Time, worldRenderView.getTime());
+			renderBlock->programParams->setFloatParameter(s_handleOcean_TileIndex, 0);
+			renderBlock->programParams->setImageViewParameter(s_handleOcean_WaveTexture, m_evolvedSpectrumTexture, 0);
+			renderBlock->programParams->setImageViewParameter(s_handleOcean_WaveTexture0, m_evolvedSpectrumTexture, 0);
+			renderBlock->programParams->endParameters(renderContext);
+
+			renderContext->compute(renderBlock);
+			renderContext->compute< render::BarrierRenderBlock >(render::Stage::Compute, render::Stage::Vertex, nullptr, 0);
+		});
+		context.getRenderGraph().addPass(rp);
+	}
+
 	for (int32_t i = 0; i < sizeof_array(m_waveTextures); ++i)
 	{
 		Ref< render::RenderPass > rp = new render::RenderPass(L"Ocean compute waves");

code/Terrain/OceanComponent.h

@@ -109,6 +109,10 @@ class T_DLLCLASS OceanComponent : public world::IEntityComponent
 	resource::Proxy< Terrain > m_terrain;
 	resource::Proxy< render::Shader > m_shaderWave;	//!< Compute shader to generate wave maps.
 	resource::Proxy< render::Shader > m_shader;
+
+	Ref< render::ITexture > m_spectrumTexture;
+	Ref< render::ITexture > m_evolvedSpectrumTexture;
+
 	Ref< render::ITexture > m_waveTextures[4];
 	Ref< const render::IVertexLayout > m_vertexLayout;
 	Ref< render::Buffer > m_indexBuffer;

data/Source/System/Shaders/Modules/Hash.xdi

@@ -4,14 +4,22 @@
 	<![CDATA[
 // Random hash functions
 
-float hash(uint n)
+float hash(uint x)
 {
-/*
-	n = (n << 13U) ^ n;
-	n = n * (n * n * 15731U + 0x789221U) + 0x1376312589U;
-	return float(n & uint(0x7fffffffU)) / float(0x7fffffff);
-*/
-	return 0.0f;
+    const uint mantissaMask = 0x007FFFFFu;
+    const uint one = 0x3F800000u;
+
+    x += (x << 10u);
+    x ^= (x >>  6u);
+    x += (x <<  3u);
+    x ^= (x >> 11u);
+    x += (x << 15u);
+
+    x &= mantissaMask;
+    x |= one;
+
+    const float r2 = uintBitsToFloat(x);
+    return r2 - 1.0f;
 }
 
 float hash12(vec2 p)

data/Source/System/Terrain/Ocean/Shaders/Wave/FFT.xdi

@@ -153,7 +153,172 @@ vec4 CalculateSpectrum(Configuration configuration, SpectrumParameters spectrum,
 	}
 }
 
+// Evolve spectrum
+
+vec2 EulerFormula(float x)
+{
+	return vec2(cos(x), sin(x));
+}
+
+vec2 ComplexMult(vec2 a, vec2 b)
+{
+	return vec2(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
+}
+
+vec4 EvolveSpectrum(Configuration configuration, ivec3 id, uint tile, vec4 initialSignal, float time, float repeatTime)
+{
+    const float lengthScale = configuration.lengthScales[tile];
+
+	//vec4 initialSignal = _InitialSpectrumTextures[uint3(id.xy, i)];
+	vec2 h0 = initialSignal.xy;
+	vec2 h0conj = initialSignal.zw;
+
+	float halfN = configuration.N / 2.0f;
+	vec2 K = (id.xy - halfN) * 2.0f * PI / lengthScale;
+	float kMag = length(K);
+	float kMagRcp = 1.0f / kMag;
+
+	if (kMag < 0.0001f)
+		kMagRcp = 1.0f;
+
+	float w_0 = 2.0f * PI / repeatTime;
+	float dispersion = floor(sqrt(configuration.gravity * kMag) / w_0) * w_0 * time;
+
+	vec2 exponent = EulerFormula(dispersion);
+	
+	vec2 htilde = ComplexMult(h0, exponent) + ComplexMult(h0conj, vec2(exponent.x, -exponent.y));
+	vec2 ih = vec2(-htilde.y, htilde.x);
+	
+	vec2 displacementX = ih * K.x * kMagRcp;
+	vec2 displacementY = htilde;
+	vec2 displacementZ = ih * K.y * kMagRcp;
+	
+	vec2 displacementX_dx = -htilde * K.x * K.x * kMagRcp;
+	vec2 displacementY_dx = ih * K.x;
+	vec2 displacementZ_dx = -htilde * K.x * K.y * kMagRcp;
+
+	vec2 displacementY_dz = ih * K.y;
+	vec2 displacementZ_dz = -htilde * K.y * K.y * kMagRcp;
+
+	vec2 htildeDisplacementX = vec2(displacementX.x - displacementZ.y, displacementX.y + displacementZ.x);
+	vec2 htildeDisplacementZ = vec2(displacementY.x - displacementZ_dx.y, displacementY.y + displacementZ_dx.x);
+        
+	vec2 htildeSlopeX = vec2(displacementY_dx.x - displacementY_dz.y, displacementY_dx.y + displacementY_dz.x);
+	vec2 htildeSlopeZ = vec2(displacementX_dx.x - displacementZ_dz.y, displacementX_dx.y + displacementZ_dz.x);
+
+	//_SpectrumTextures[uint3(id.xy, i * 2)] = vec4(htildeDisplacementX, htildeDisplacementZ);
+	//_SpectrumTextures[uint3(id.xy, i * 2 + 1)] = vec4(htildeSlopeX, htildeSlopeZ);
+
+	return vec4(htildeDisplacementX, htildeDisplacementZ);
+}
+
+// FFT
+
+vec2 ComplexExp(vec2 a)
+{
+    return vec2(cos(a.y), sin(a.y) * exp(a.x));
+}
+
+vec4 ComputeTwiddleFactorAndInputIndices(Configuration configuration, uvec2 id)
+{
+    uint b = configuration.N >> (id.x + 1);
+    vec2 mult = 2 * PI * vec2(0.0f, 1.0f) / configuration.N;
+    uint i = (2 * b * (id.y / b) + id.y % b) % configuration.N;
+    vec2 twiddle = ComplexExp(-mult * ((id.y / b) * b));
+    return vec4(twiddle, i, i + b);
+}
+
+#define SIZE 512
+#define LOG_SIZE 9
+
+shared vec4 fftGroupBuffer[2][SIZE];
+
+void ButterflyValues(uint step, uint index, out uvec2 indices, out vec2 twiddle)
+{
+    const float twoPi = 6.28318530718;
+    uint b = SIZE >> (step + 1);
+    uint w = b * (index / b);
+    uint i = (w + index) % SIZE;
+
+    // sincos(-twoPi / SIZE * w, twiddle.y, twiddle.x);
+	twiddle.x = cos(-twoPi / SIZE * w);
+	twiddle.y = sin(-twoPi / SIZE * w);
+
+    // This is what makes it the inverse FFT
+    twiddle.y = -twiddle.y;
+    indices = uvec2(i, i + b);
+}
+
+vec4 FFT(Configuration configuration, uint threadIndex, vec4 inputData)
+{
+    fftGroupBuffer[0][threadIndex] = inputData;
+    
+	// GroupMemoryBarrierWithGroupSync();
+	groupMemoryBarrier();
+	barrier();
+
+    int flag = 0;
+
+    [[unroll]]
+    for (uint step = 0; step < LOG_SIZE; ++step)
+	{
+        uvec2 inputsIndices;
+        vec2 twiddle;
+        ButterflyValues(step, threadIndex, inputsIndices, twiddle);
+
+        vec4 v = fftGroupBuffer[flag][inputsIndices.y];
+        fftGroupBuffer[1 - flag][threadIndex] = fftGroupBuffer[flag][inputsIndices.x] + vec4(ComplexMult(twiddle, v.xy), ComplexMult(twiddle, v.zw));
+
+        flag = 1 - flag;
+
+        // GroupMemoryBarrierWithGroupSync();
+		groupMemoryBarrier();
+		barrier();
+    }
+
+    return fftGroupBuffer[flag][threadIndex];
+}
+
 // Displacement
+/*
+void GenerateDisplacement(Configuration configuration, ivec3 id, uint tile)
+{
+	const int i = 0;	// Tile
+
+    vec4 htildeDisplacement = Permute(_SpectrumTextures[uvec3(id.xy, i * 2)], id);
+    vec4 htildeSlope = Permute(_SpectrumTextures[uvec3(id.xy, i * 2 + 1)], id);
+
+        vec2 dxdz = htildeDisplacement.rg;
+        vec2 dydxz = htildeDisplacement.ba;
+        vec2 dyxdyz = htildeSlope.rg;
+        vec2 dxxdzz = htildeSlope.ba;
+        
+        float jacobian = (1.0f + _Lambda.x * dxxdzz.x) * (1.0f + _Lambda.y * dxxdzz.y) - _Lambda.x * _Lambda.y * dydxz.y * dydxz.y;
+
+        vec3 displacement = float3(_Lambda.x * dxdz.x, dydxz.x, _Lambda.y * dxdz.y);
+
+        vec2 slopes = dyxdyz.xy / (1 + abs(dxxdzz * _Lambda));
+        float covariance = slopes.x * slopes.y;
+
+        float foam = _DisplacementTextures[uint3(id.xy, i)].a;
+        foam *= exp(-_FoamDecayRate);
+        foam = saturate(foam);
+
+        float biasedJacobian = max(0.0f, -(jacobian - _FoamBias));
+
+        if (biasedJacobian > _FoamThreshold)
+            foam += _FoamAdd * biasedJacobian;
+
+        _DisplacementTextures[uint3(id.xy, i)] = float4(displacement, foam);
+        _SlopeTextures[uint3(id.xy, i)] = float2(slopes);
+
+//        if (i == 0) {
+//            _BuoyancyData[id.xy] = displacement.y;
+//        }
+    }
+}
+*/
+// Old
 
 float CalculateDisplacement(vec2 position, float time)
 {
@@ -168,7 +333,6 @@ float CalculateDisplacement(vec2 position, float time)
 	return v;
 }
 
-
 	]]>
 	</text>
 </object>