Dear Gamedev.net-aratti

I've spent the last week or so rendering a simple ocean using gerstner waves but having issues with tiling, so I decided to start rendering them "properly" and dip my toes into the murky waters of rendering a heightfield using an iFFT.

There are plenty of papers explaining the basic gist -

1) calculate a frequency spectrum

2) use this to create a heightfield using ifft to convert from frequency domain to spatial domain - animating with time t

Since the beginning of this journey I have learned about things like the complex plane, the complex exponent equation, the FFT in more detail etc but after the initial steps of creating an initial spectrum (rendering a texture full of guassian numbers with mean 0 and sd of 1, filtered by the phillips spectrum) I am still totally lost.

My code for creating the initial data is here (GLSL):

float PhillipsSpectrum(vec2 k){
 
  //kLen is the length of the vector from the centre of the tex
  float kLen = length(k);
  float kSq = kLen * kLen;
 
  // amp is wave amplitude, passed in as a uniform
  float Amp = amplitude;
 
  //L = velocity * velcoity / gravity
  float L = (velocity*velocity)/9.81;
 
  float dir = dot(normalize(waveDir),normalize(k));
 
  return Amp * (dir*dir) * exp(-1.0/(kSq * L * L))/ (kSq * kSq) ;
}
 
void main(){
 
   vec3 sums;
 
   //get screenpos - center is 0.0 and ranges from -0.5 to 0.5 in both
   //directions
   vec2 screenPos = vec2(gl_FragCoord.x,gl_FragCoord.y)/texSize - vec2(0.5,0.5);
 
   //get random Guass number
   vec2 randomGuass = vec2(rand(screenPos),rand(screenPos.yx));
 
   //use phillips spectrum as a filter depending on position in freq domain
   float Phil = sqrt(PhillipsSpectrum(screenPos));
   float coeff = 1.0/sqrt(2.0);
 
   color = vec3(coeff *randomGuass.x * Phil,coeff * randomGuass.y * Phil,0.0);
}

which creates a texture like this:

[attachment=26468:wave.png]

Now I am totally lost as how to :

a) derive spectrums in three directions from the initial texture

b) animate this according to time t like mentioned in this paper (https://developer.nvidia.com/sites/default/files/akamai/gamedev/files/sdk/11/OceanCS_Slides.pdf) on slide 5

I might be completely stupid and overlooking something really obvious - I've looked at a bunch of papers and just get lost in formulae even after acquainting myself with their meaning. Please help.

From reading Tessendorf's paper and all the other code examples out there (and doing an implementation in Matlab), I think I can try to offer some explanations.

It seems you're getting the correct \(H_0\) spectrum from the texture. My image (in grayscale) and images in other papers look similar to yours:

[attachment=26503:philspect.PNG]

To get the height map, you need to use your texture values in the following formula and you can plug in the time parameter there as well:

\[H(k,t) = H_0(k,t)e^{i \omega t} + H_0^* e^{-i \omega t},\, \omega = g \left | k \right |\]

where \(H_0^*\) is the complex conjugate of \(H_0\). I think since your Y coordinate in the color vec3 seems to be the "imaginary" component of the complex number, then you can get the complex conjugate by just multiplying Y by -1. This \(H(k,t)\) is your height map. Here you can plug in your time \(t\) and get the animated displacements (at least I think it works; I didn't do the full animation). The decomposition into the spatial directions is given by:

\[D_x(k)=i\frac{k_x}{\left | k \right |}H(k),\,D_y(k)=i\frac{k_y}{\left | k \right |}H(k)\]

From there, you take the inverse FFT to get your spatial coordinates. In my Matlab implementation, those came out as complex numbers so I just took the norm and used that for the values, but I think you can just use the real parts as well. My 3 spectrum images look like this (dx, height, dy):

[attachment=26504:philspect2.PNG]

I haven't gotten around to using these to actually create waves. I've never done water simulations before, but the problem interested me from a theoretical standpoint so I figured I'd try it out. I'll post my Matlab code and if you use NumPy/SciPy the functions I used are very similar if you want to load it up and mess with it. Hope this helps!

clear;
clc;
close all;

M = 512;    % resolution of grid
N = 512;    % resolution of grid
L = 2000;   % height of waves used to calculate dx and dz

[m,n] = meshgrid(-M/2:1:M/2, -N/2:1:N/2);
kx = 2*pi/L * m;
ky = 2*pi/L * n;

W = [30,0];        % horizontal wind vector (m/s)
V = norm(W);        % m/s?
What = W / V;
g = 9.81;           % m/s^2
A = 1.0;            % amplitude of wave

% compute Phillips spectrum and freq amplitudes
P = zeros(M,N);
H0 = zeros(M,N);
for i=1:M
    for j=1:N
        k = [kx(i,j),ky(i,j)];
        nk = norm(k);
        if(nk == 0)
            mag = 0;
        else
            khat = k / nk;
            mag = A / nk^4 * exp(-1/(nk * V^2/g)^2) * ((dot(khat,What)))^2;
        end
        damping = exp(-nk^2*(L/1000)^2);
        val = mag * damping;
        P(i,j) = val; 
        H0(i,j) = 1/sqrt(2) * (randn + randn * 1i) * sqrt(val);
    end
end

figure(1)
    imshow(abs(H0))

% compute heightmap
H = zeros(M,N);
Dx = zeros(M,N);
Dy = zeros(M,N);
t = 1;
for i=1:M
    for j=1:N
        k = [kx(i,j),ky(i,j)];
        nk = norm(k);
        w = sqrt(g*nk);
        val = H0(i,j)*exp(1i*w*t) + conj(H0(M-i+1, N-j+1))*exp(-1i*w*t);
        H(i,j) = val;
        if(nk == 0)
            Dx(i,j) = 0;
            Dy(i,j) = 0;
        else
            Dx(i,j) = 1i*k(1)/nk * val;
            Dy(i,j) = 1i*k(2)/nk * val;
        end
    end
end

dx = abs(ifft2(Dx));
dy = abs(ifft2(Dy));
dz = abs(ifft2(H));

figure(2)
    subplot(1,3,1)
    imshow(dx)
    subplot(1,3,2)
    imshow(dz)
    subplot(1,3,3)
    imshow(dy)

Thank you very much. I saw the Keith Lantz links before but all the derivations and mathematical breakdowns made my head hurt. I can understand your explanations a lot better! i.e "complex conjugate in this case is multiplying by -1". If I have any more problems I'll let you know! Thanks again