# FFT Water Video shows strange behaviour

This topic is 1463 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

Hi all,

I've implemented radix-2 FFT on the gpu and that's the result that I've got so far (for the lod algorithm I'm using geometry clipmap with toroidal update on).

sometime It looks like that the waves instead of doing a smooth and gentle fade off, they just looks like they cut off with a discontinuity ... If you look at the 1st video at 0:18 I show the heightmap update along with the related spectrum. The heightmap looks ok to me (it actually has choppiness in as well, as you can see because it's colored).

Note: there is choppiness as well (which I tried to disable it to see if that was the problem, but it wasn't).

can someone shed some light or has some ideas ?

Also what could be a good rule of thumb to follow to try different settings for the water ?

Reminder:

I've got grid size in meters (MxN) and grid resolution which is the actual heightmap res.

wind speed V

wind direction W

PS.: if I increase V the speed of the simulation is still the same. I guess V will generate just bigger waves ? ...

##### Share on other sites

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.

##### Share on other sites

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  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
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);
float2 Knorm         = normalize(K);

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

//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 ?

Edited by MegaPixel

##### Share on other sites

R16G16B16A16_FLOAT is an appropriate format although I would use an integer format as it provides a more uniform distribution of values. For that you'd need to know ahead the bounds of your output data.

Your calculation of omega appears correct. The lateral displacement is wrong instead. Here's the correct formula:

float2 DxKt = HKt.yx * float2(1, -1) *  K.x / K = HKt.yx * float2(1, -1) * Knorm.x;

float2 DyKt = HKt.yx * float2(1, -1) *  K.y / K = HKt.yx *float2(1, -1) * Knorm.y;

Basically you flipped the sign of the lateral components. Think about Gerstner waves, characterized by a 90 degrees phase offset between the vertical and lateral displacement. For a complex wave that translates to the transformation from  (cos + i sin) to (sin - i cos)

May I be critical and ask you about your goal ? Copying and pasting code from other sources is dull and pointless. In my opinion it is essential to have a solid grasp of the underlying theory before moving to an implementation, especially on advanced topics like this.

##### Share on other sites

R16G16B16A16_FLOAT is an appropriate format although I would use an integer format as it provides a more uniform distribution of values. For that you'd need to know ahead the bounds of your output data.

Your calculation of omega appears correct. The lateral displacement is wrong instead. Here's the correct formula:

float2 DxKt = HKt.yx * float2(1, -1) *  K.x / K = HKt.yx * float2(1, -1) * Knorm.x;

float2 DyKt = HKt.yx * float2(1, -1) *  K.y / K = HKt.yx *float2(1, -1) * Knorm.y;

Basically you flipped the sign of the lateral components. Think about Gerstner waves, characterized by a 90 degrees phase offset between the vertical and lateral displacement. For a complex wave that translates to the transformation from  (cos + i sin) to (sin - i cos)

May I be critical and ask you about your goal ? Copying and pasting code from other sources is dull and pointless. In my opinion it is essential to have a solid grasp of the underlying theory before moving to an implementation, especially on advanced topics like this.

I don't understand why you flipped the immaginary part and real part of H(K,t) ... (in the paper is not clear how that calculation is done ... and it's just says -i * Knorm * H(K,t), to me looks like it's that -i is a complex number with no real part and just Im(-i) = -1 and I can't see any flip in the immaginary and real part as you did ...). I'm sure I'm missing something from the equation that I was trying to understand

I didn't copy and paste the code (the pixel shaders that you see there are all written by myself from scratch straight from theory), I actually studied quite a lot before to implement it and then implemented everything by myself from scratch . There is not a special goal, it's

just for learning purposes ;-).

I'll try the correction and thanks for the explanation ;). Also, your water is beautiful

Edited by MegaPixel

##### Share on other sites

R16G16B16A16_FLOAT is an appropriate format although I would use an integer format as it provides a more uniform distribution of values. For that you'd need to know ahead the bounds of your output data.

Your calculation of omega appears correct. The lateral displacement is wrong instead. Here's the correct formula:

float2 DxKt = HKt.yx * float2(1, -1) *  K.x / K = HKt.yx * float2(1, -1) * Knorm.x;

float2 DyKt = HKt.yx * float2(1, -1) *  K.y / K = HKt.yx *float2(1, -1) * Knorm.y;

Basically you flipped the sign of the lateral components. Think about Gerstner waves, characterized by a 90 degrees phase offset between the vertical and lateral displacement. For a complex wave that translates to the transformation from  (cos + i sin) to (sin - i cos)

May I be critical and ask you about your goal ? Copying and pasting code from other sources is dull and pointless. In my opinion it is essential to have a solid grasp of the underlying theory before moving to an implementation, especially on advanced topics like this.

I don't understand why you flipped the immaginary part and real part of H(K,t) ... (in the paper is not clear how that calculation is done ... and it's just says -i * Knorm * H(K,t), to me looks like it's that -i is a complex number with no real part and just Im(-i) = -1 and I can't see any flip in the immaginary and real part as you did ...). I'm sure I'm missing something from the equation that I was trying to understand

I didn't copy and paste the code (the pixel shaders that you see there are all written by myself from scratch straight from theory), I actually studied quite a lot before to implement it and then implemented everything by myself from scratch . There is not a special goal, it's

just for learning purposes ;-).

I'll try the correction and thanks for the explanation ;). Also, your water is beautiful

I verified the calculation, now I understand why was uncorrect. I forgot to consider that there will be i^2 after the multiplication that obviously is -1 :) and the immaginary part will appear on the real part which then is meant to be swapped again because I store the real part in the first component and the immaginary part in the second ;)

##### Share on other sites

I'm glad to have been of help. In my experience it helps immensely to understand the physical meaning of all terms of a mathematical model. Sometimes a simple change of convention, untold assumptions or misprints in a publication lead to bugs that are very hard to fix when you are unfamiliar with the model.

Was the bizarre animation of your waves caused by the flipped sign in the lateral displacement ? I would like to see an updated video of your technique and if possible a demo.

##### Share on other sites

I'm glad to have been of help. In my experience it helps immensely to understand the physical meaning of all terms of a mathematical model. Sometimes a simple change of convention, untold assumptions or misprints in a publication lead to bugs that are very hard to fix when you are unfamiliar with the model.

Was the bizarre animation of your waves caused by the flipped sign in the lateral displacement ? I would like to see an updated video of your technique and if possible a demo.

Yes that's true. But is also true that if the main sources of the mathematical model are no very precise and clear in stating what assumptions are made and so on, someone

that has to start studying them might not be find himself in the right track from the beginning. Also, the two main sources that I've checked are obviously the J. Tessendorf paper and the gamasutra article (which is actually a summarization of approaches) from Jensen. I understand that both paper have advanced mathematical models, like the DFT to work in the frequency domain and back to time domain etc.

If you know about some other papers that are less advanced that I can study to get an introduction (other than what I've studied in my engineering degree) I'll be happy to study them as well in order to integrate my knowledge.

Btw, there was a bug in the displacement of the y axis (it was flipped as well). Also the horizontal disp :). I'll have a look again at the papers to see if I've missed something.

In the meantime thanks a lot for your advices and help and keep up with your awesome work.

Updated video: