What texture format do you use for storing the generated heightmap ? If it's an integer format, then it is possible that the results are clamped. Also, the wave shape looks off somehow, as if the lateral displacement was not in sync with the vertical one. Some code would help here to better understand the issue.

The wind speed affects the waves amplitude based on the spectrum you are using (Philips, JONSWAP, etc ?). The speed of a wave depends only on its spatial frequency, according to the wave dispersion relation.

Dx, Dy, Dz in the space domain (after inverse FFT has been applied) are stored in R16_FLOAT and then interleaved in a R16G16B16A16_FLOAT for ease of use later on (so no integer here but I suppose that 16bits float is enough here ? o should I use 32bit ?.

I use Phillips spectrum for now.

The speed of wave dependes on the dispersion relation ? you mean omega? at the moment I'm using omega = sqrt(k*g) trying to let it loop with period T.

For the choppiness displacement vector I'm not sure if it's correctly calculated, because on the tessendorf paper as well as the nvidia slides is not very clear how to carry on

the calculation. Also shouldn't the horizontal displacement vector be the gradient vector flipped ?

this is the pixel shader code to generate H0 before to start the simulation (this will be execute at the simulation start and any time that a parameter like wind speed, wind direction etc. will change).

float lengthSqr(in float2 v){
return dot(v,v);
}
float getPhillipsSpectrum(in float2 k,in OceanData oceanData){
const float windDependency = oceanData.A_V_windDependency_unused.z;
const float2 W = oceanData.W_unused_unused.xy;
const float A = oceanData.A_V_windDependency_unused.x;
const float V = oceanData.A_V_windDependency_unused.y;
float L = (V*V)/GRAVITY;
float damping = 0.001f;
float l = L*damping;
float ksqr = lengthSqr(k);
if(ksqr < 0.0000001f) //avoid division by 0
return 0.f;
float2 Wn = normalize(W);//normalize wind direction
float2 kn = normalize(k);
float kdotw = dot(kn,Wn);
float k4 = ksqr*ksqr;
float kL2 = ksqr*L*L;
float exp_term = exp(-1.f/kL2);
float P_k = A*(exp_term/(ksqr*ksqr))*(kdotw*kdotw); //resulting Phillips spectrum
//introduce wind dependency
if(kdotw < 0.00001f)
P_k *= windDependency;
//finally suppress waves smaller than a small length (l<<L)
return P_k*exp(-ksqr*l*l);
}
float2 get_k(in float2 pos,in OceanData oceanData){
float2 k;
const float4 width_height_Lx_Lz = oceanData.width_height_Lx_Lz;
k.y = ((pos.y-width_height_Lx_Lz.y*0.5f)*TWOPI)/width_height_Lx_Lz.w;
k.x = ((pos.x-width_height_Lx_Lz.x*0.5f)*TWOPI)/width_height_Lx_Lz.z;
return k;
}
float4 main(PS_Input Input) : SV_Target{
float2 pos = Input.Pos.xy; //screen space position
float2 gaussRand0 = gaussianDistributionTexture.Load(int3(pos,0)).xy;
float2 gaussRand1 = gaussianDistributionTexture.Load(int3(pos,0)).zw;
float2 k = get_k(pos,oceanData);
//get the Phillips spectrum
float P_k = sqrt(getPhillipsSpectrum(k,oceanData)*0.5f);
float P_mk = sqrt(getPhillipsSpectrum(-k,oceanData)*0.5f);
//H0(k)
float2 H0k = P_k*gaussRand0;
//H0*(-k)
float2 H0mkConj = float2(P_mk*gaussRand1.x,-P_mk*gaussRand1.y);
return float4(H0k,H0mkConj);
}

this is my pixel shader code for the three spectrum update:

PS_Output main(PS_Input Input){
PS_Output Output = (PS_Output)0;
float2 pos = Input.Pos.xy; //screen space position
float4 h0k_h0mkConj = h0k_h0mkConjTexture.Load(int3(pos,0));
float2 K = get_k(pos,oceanData); //get the wave vector
//start calculating H(k,t) first then Dx and Dz
const float T = oceanData.A_V_windDependency_T.w;
const float t = time.x;
const float2 H0k = h0k_h0mkConj.xy;
const float2 H0mkConj = h0k_h0mkConj.zw;
float omega0 = TWOPI/T; //basic angular frequency
float k = length(K);
float omegad = floor(sqrt(k*GRAVITY)/omega0)*omega0;
float2 Knorm = normalize(K);
float omegadt = omegad*t;
float2 expLeft = float2(cos(omegadt),sin(omegadt)); //euler formula
float2 expLeftConj = conjComplex(expLeft);
//H(k,t)
//HKt = H0*expLeft + H0Conj*conj(expLeft);
float2 a = mulComplex(H0k,expLeft);
float2 b = mulComplex(H0mkConj,expLeftConj);
//Dy axis displacement in frequency domain
float2 HKt = addComplex(a,b);
//Dx,Dz displacements in frequency domain
float2 DxKt = mulComplex(float2(0.f,-Knorm.x),HKt);
float2 DzKt = mulComplex(float2(0.f,-Knorm.y),HKt);
Output.Dx = float4(DxKt,0.f,1.f);
Output.Dy = float4(HKt,0.f,1.f);
Output.Dz = float4(DzKt,0.f,1.f);
return Output;
}

those two are the complex representation of the horizontal displacement vector in frequency domain:

float2 DxKt = mulComplex(float2(0.f,-Knorm.x),HKt); //mulComplex is used to perform the complex multiplication between two complex number

float2 DzKt = mulComplex(float2(0.f,-Knorm.y),HKt);

I'm not sure if the computation is correct as in the paper it doesn't break down further so I would expect to be like (the nvidia slides don't use the minus sign ...):

-i * (K/k) * H(K,t)

which means: -i is a complex number with no real part but just Im(-i) = -1 if I understand correctly.

so this should come down to:

Kn == K/k which is K normalized

-i*Kn.x ---> float2(0.f,-Knorm.x)

-i*Kn.y ---> float2(0.f,-Knorm.y)

then to produce each component I should perform a complex multiplication with H(k,t) for each component to produce Dx and Dz in the frequency domain.

I paste you the nvidia approach (which is actually the tessendorf one):

Also: nvidia doesn't put any minus sign while Tessendorf does ... !!! ... :/

Plus, the colors of the final displacement from NVIDIA seem to be differents in general from the ones of my heightmap as you can see in my first video ... might be because of some bug in the horizontal vec displacement ?

thanks in advance for your help

**Edited by MegaPixel, 15 January 2014 - 11:36 AM.**