Opengl Texture LoD

LoD selection by opengl

When rendering a textured object, opengl applies filtering on the texture to avoid aliasing artefacts. This is why the texture is loaded in a Mipmap pile. The base level of this pile is the original texture and then every level is built by averaging neighbours texels to create a texture of half the size of the previous one.


Then, during rendering, opengl picks the texture from the right LoD. The choice of this LoD is governed by a scale factor  scale_factor. This factor is used to select the LoD level λ






The bias parameters are either defined by the TEXTURE_LOD_BIAS parameter or by the fragment shader. In λ’, those bias are clamped between -MAX_TEXTURE_LOD_BIAS and MAX_TEXTURE_LOD_BIAS.

This leaves the question of the scale factor scale_factor. It is defined as


Illustrated in 2D


Where P(u,v, w) are projected texture coordinates and p(x, y) image coordinates. When using a 1D or 2D texture, the corresponding derivatives are set to 0.

This relies on the fact that when the object is far away from the camera, the point projecting on a pixel will be distant for the points projecting on the neighbouring pixels.

Biasing a LoD

There is therefore two ways to bias a LoD. Either by changing the gradients, therefore altering scale_factor, or by adding a fixed bias.

Fixed Bias

Here, setting the bias during a texture call in the fragment shader or by setting the TEXTURE_LOD_BIAS parameter makes no difference. In our example, we used

texture(sampler, coord, [bias])

in the fragment shader.

Without any bias, we get

color = texture(tex, v_texcoord, 0);
Texture from


If we set a bias to -4, every lod level chosen will be offset by -4. Therefore, when opengl computes lod = 0, 1, 2, 3 or 4, it will use the level lod – 4 = 0.

color = texture(tex, v_texcoord, -4);
Texture from


Oppositely, if we set a bias to 4, the lod levels will be offset by 4, leading to coarser LoD being used.

color = texture(tex, v_texcoord, 4);
Texture from


Changing the gradients

This can be done using the

textureGrad(sampler, coord, dPdx, dPdy)

where dPdx and dPdy and the derivative of the texture coordinates when looking at the neighbouring pixel respectively on the x and y axis.

For a use case, see the Screen Projected Grid post.

Enforcing a LoD

You can force the LoD selection using the function

textureLod(sampler, texture coordinate, lod)

When using the simple texture function, the correct LoD level is chosen automatically by opengl. For example, when rendering an infinite plane, the lod used in the background is not the same as in the foreground.

Texture from


If we render the same scene using

color = textureLod(tex, v_texcoord, 0);

to force opengl to always use the base texture level, we get

Texture from


We clearly notice aliasing in the background. If we do it the other way around, rendering it using always a high Lod level, we get

color = textureLod(tex, v_texcoord, 6);
Texture from

Since we are using a texture that is a far average of the original one, the resulting textured plane is heavily blurred.


[1] Opengl Specification Version 4.5 –

Screen Space Grid

Rendering terrains or ocean is usually done by displacing vertices of a grid defined in world space. Depending on the resolution of the grid, this might create aliasing artefacts when several triangles project onto the same pixel and can be computationally expensive. This is why Level of Details algorithms are used. Those methods consist in adapting the resolution of the grid to its distance to the camera.

Example of LoD grid.


Another solution would be to define the grid in screen space and then to back project it onto an infinite plane. With such techniques, the grid resolution on screen in constant and the change in resolution in world space is continuous (no hard transition between different resolutions).

Grid defined in screen space and back projected on a plane. Image from [1]

This algorithm works by

  • Creating a grid orthogonal to the viewer in screen space.
  • Project it on a plan in world space
  • Displace the points according to a height field.
  • Render the projected plane

In its naive form, this algorithm fails for each point of the far plane of the frustum that does not intersect the plane.

When the camera looks above the horizon, part of the frustum points will never intersect the plane: the intersection points will computed behind the camera. Image from [2]
Furthermore, it is important that the camera is sufficiently high upon the plane to prevent the points after displacement to appear below it

visible_rangeThose are heavy constraints on the camera moves.

First, we create a grid in screen space. Therefore, we generate a grid of vertices with 3D coordinates ranging from [-1, -1, 0] to [1, 1, 0]. If you visualize them with

#version 330


layout (location=0) in vec3 position;

void main()
    gl_Position = vec4( position.x, position.y, 0, 1 );



uniform sampler2D tex;
uniform int line;

out vec4 color;

in vec3 test;

void main()
    color = vec4(0, 0, 0, 1);


You get something like this (in wireframe)


Then, each of those points has to be projected onto a 2D plane.


Here, we want to find the point P, projection of p on the world plane. To do so, we will use a simple Pythagore. The triangle made by p, P and h is a squared triangle. the distance between h and p is the height of the camera. The vector from p to P is the world direction coming from the camera through the vertex p. Then, to  get P, all we need is the distance t and to solve

P = p + world direction * t

Pythagore tells us that

t = cos(α) / camera_height

with cos(α) = dot(up_vector, world_dir).

This leads to the code ( my up vector is vec3(0, -1, 0) ).

vec3 toWorldPos(vec3 posScreen)
    vec4 vertex = vec4(posScreen, 1);
    vec3 camera_dir = 
          normalize( ( inverse(projection) * vertex ).xyz );
    vec3 world_dir = 
          normalize( inverse(view) * vec4(camera_dir, 0) ).xyz;

    float t = camera_position.y / -world_dir.y;

    return camera_position + t*world_dir;

And now , you can displace the computed world position with a heightmap and render those vertices with the usual Model View Projection transforms. Here, we use as a heightmap


#version 330


uniform mat4 projection;
uniform mat4 view;
uniform vec3 camera_position;

uniform sampler2D tex;

layout (location=0) in vec3 position;

out vec3 world_pos;

void main()
    world_pos = toWordPos( vec3(position.x, 0, position.y) );
    float height = texture2D(tex, world_pos.xz).r;
    test.y = height;
    gl_Position = projection * view * vec4(world_pos, 1);




out vec4 color;

in vec3 world_pos;

void main()
    vec3 normal = normalize(cross( 
    vec3 light = vec3(0, 1, 1);
    vec3 k = vec3(1, 1, 1);
    color = vec4(1, 0, 0, 0) * 
            dot(normalize(-light), normalize(normal));


Leading to



huuuum this wireframe is a bit fishy ain’t it ?

Actually, this method looks nice and easy but it has a major flaw; when there is no intersection between the world direction and the plane is completely undefined. Therefore, if you have a camera looking along the z axis, half the screen is not intersecting with the plane.

A solution to this problem is to project the points in a geometry shader. This way, we will be able to discard the triangles that are not in front of the camera.


layout (location=0) in vec3 position;

out vec3 vertex_position;

void main()
    vertex_position = position;
    gl_Position = vec4(position, 1);




layout(triangles) in;
layout(triangle_strip, max_vertices = 3) out;

uniform mat4 projection;
uniform mat4 view;
uniform vec3 camera_position;

uniform sampler2D tex;

in vec3 vertex_position[];

out vec3 geometry_position;

bool toWorldPos(vec3 posScreen, out vec3 pos)
    vec4 vertex = vec4(posScreen, 1);
    vec3 camera_dir = 
             normalize( ( inverse(projection) * vertex ).xyz );
    vec3 world_dir = 
             (inverse(view) * vec4(camera_dir, 0)).xyz;
    if (world_dir.y == 0)
    return false;

    float t = camera_position.y / -world_dir.y;
    if (t < 0)
    return false;
    pos = camera_position + t*world_dir;
    pos.y = 0;
    return true;

void main() {
    float fact = 0.01;
    for(int i=0; i<3; i++)
      if(toWorldPos(vertex_position[i], geometry_position))
        vec4 disp = texture(tex, geometry_position.xz);
        geometry_position = vec3(geometry_position.x, 
        gl_Position = 
              projection * view * vec4(geometry_position, 1);

Leading to this result


Until there everything seems fine … then why is the plane so ugly ???

Actually it is once again a matter of filtering. For opengl, the grid is right in front of the camera. Therefore, it always accesses the texture to its highest level of detail. Quickly, the grid resolution is not enough to sample the texture and aliasing occur.

This is the same object, rendered using a 50×50 and a 150×140 screen grid.


If you want to access the correct LoD, you need to use the function

textureGrad(samples, coord, gradX, gradY)

This function changes the way opengl selects its LoD.For detailed information about opengl LoD selection process, see the post OpenGL Texture Access.

Here, it requires computing the derivative between the projection of one point of the projective grid and the projections of its neighbours on X and Y.

vec3 position_x;
vec3 position_y;
toWorldPos( vertex_position[i] + vec3(1.0/150.0, 0, 0),  
toWorldPos(vertex_position[i] + vec3(0, 1.0/150.0, 0), 
vec2 dudx = position_x.xz - geometry_position.xz;
vec2 dudy = position_y.xz - geometry_position.xz;
//We multiply by 0.1 to have less texture tiles
vec4 disp = textureGrad(tex, 
                       dudx*0.1, dudy*0.1);

Which gives us as a final rendering


[1] J. D. Mahsman, Projective grid mapping for planetary terrain. University of Nevada, Reno, 2010.
[2] “Real-time water rendering – Introducing the projected grid concept – Habib’s Water Shader.”

Gooch Shading

Photorealistic rendering is a major issue nowadays. We are looking for increased realism in moves and games. But this realism can come as a flaw when rendering for illustration purposes. This lead to a whole field of study called NPR, Non Photorealistic Rendering. An important work in this field is the shading of Gooch [1]. This is a method to shade objects that focuses on rendering an picture that is easily readable.

The first point made by Gooch is that the edges of the objects should be drawn. There are several methods to do so. The fastest one to implement is to use two rendering passes. The first one draw the wireframe of the object with very large lines. The second one draws the object itself. The depth test hides the edges that are not visible. This methods is easy to implement but the results present some artefacts where edges connect.

The edges help to determine the silhouette of the skeleton. However, you may note that since the edges are drawn in black, it can be difficult to discern the edges from the shadowed part of the mesh. T rex mesh from :


You can note that the lines are broken where edges connect. Those artefacts could be avoided by using another method for the edge detection and silhouette drawing.


glCullFace (GL_FRONT);
glDepthFunc (GL_LEQUAL);
glPolygonMode (GL_FRONT_AND_BACK, GL_LINE); 
glCullFace (GL_BACK);
glDepthFunc (GL_LESS);

Once the edges are drawn, the second point made by Gooch is that the shadows and highlights on the object should be of different color than the edges or the object. Otherwise, different parts of the object can be rendered in the same color which prevents them from being easily differentiable.

Therefore, what is suggested is to take a « hot » color (like orange) for the enlightened parts of the object, and a « cold » color (like blue) for the shadowed parts. This way, the object colors interpolate between orange and blue, and therefore are never white or black. The edges and the highlights of the object are then perfectly readable.

In his paper, Gooch incorporates the original color of the object. The « hot » and « cold » colors thus become


where kd is the color of the object, kblue is (0, 0, 1) and kyellow is (1, 1, 0), α = 0.2 and β = 0.6.

Leading to the final shading formula


where l is the normalized light direction and n is the normalized normal.

There is the same mesh enlightened with Phong and with Gooch.

T Rex lighten with Phong. T rex mesh from :
The same T Rex with gooch shading. It is much easier to discern the different bones. T rex mesh from :
vec4 gooch_shading(vec4 m_color, //color of the mesh
                   float m_shine, //shininess of the surface
                   vec3 l_direction, //light direction
                   vec3 v_normal, //normal
                   vec3 c_direction) //camera direction
    float kd = 1;
    float a = 0.2;
    float b = 0.6;

    float NL = dot(normalize(v_normal), normalize(l_direction));
    float it = ((1 + NL) / 2);
    vec3 color = (1-it) * (vec3(0, 0, 0.4) + a* 
               +  it * (vec3(0.4, 0.4, 0) + b*;
    vec3 R = reflect( -normalize(l_direction), 
                      normalize(v_normal) );
    float ER = clamp( dot( normalize(c_direction), 
                     0, 1);
    vec4 spec = vec4(1) * pow(ER, m_shine);

    return vec4(, m_color.a);
[1] A.Gooch, B. Gooch, P. Shirley, and E. Cohen, “A non-photorealistic lighting model for automatic technical illustration,” in Proceedings of the 25th annual conference on Computer graphics and interactive techniques, 1998, pp. 447–452.

Ocean modelisation – Fourier

Modeling ocean waves can be a difficult task. Waves are generated from the interaction between the wind and the ocean surface and this interaction is usually overly simplified during modelisation. There are two main classes of solution to generate waves, either using the Fourier spectra of the ocean [2], or representing the waves as a sum of trochoids [1].

Modeling terrains using Fourier is quite known. You get the Fouier Transform of a white noise (so that all frequencies appear), you filter it to select a subset of those frequencies and you use the inverse Fourier Transform on this subset. This gives you a heightmap representing your terrain ( — Article on « Frequency synthesis of Landscapes (and clouds) » ).

The same process can be applied to the oceans. When modeling the ocean using a Fourier spectra, the idea is to represent a set of waves as a sum of sinusoids with different frequencies. The trick is therefore to choose those frequencies carefully. Oceanographic research have actually an answer to this as several spectrum where defined to represent the frequency of waves. By applying a Fourier transform on a subset of those frequency, you can get an heightmap representing a very realistic ocean.

Here, we will used the formulas given by [3]. Those formulas and the corresponding code are already very well explained in the tutorial by Keith Lantz ( but i’ll reexplain it anyway.

The idea is to disturb a grid mesh of M*N points.


Those points will be displaced to form the ocean surface. The question is how to displace them. It all starts with an inverse Fourier Transform


where h(x, t) is the height for a point x of coordinates (x, y, z) at a time t. k are the frequencies of the complex exponentials and h~ is the Fourier Transform of h. Since we need as many frequencies as the number of points, we put k = ( (2πn / Lx ), (2πm / Lz) ) with n in [-N/2, N/2] and m in [-M/2, M/2].

To compute h with a discrete fourier transform, we need to sum the computed h~. We have


with ω(k)² = gk, g being the gravitational constant and k being the length of the vector k. This leaves the question of evaluating the complex number h0~, h0~* being its conjugate. Expressing h~ as a sum allows to maintain the complex conjugation property, h~(k, t) = h~*(-k, t).

We know that the waves follow a spatial spectrum similar as the Phillips spectrum for wind driven waves


where L = V²/g; V being the windspeed and w being the direction of the wind. The ^ symbol means that the vectors are normalized. A is a constant (0.0001 in our examples).

This formula is furthermore multiplied by


with l = 0.001 * L. This removes the waves too small to be seen.

From there, our h0~ takes samples from this spectrum as


where ξ are independent draws from a gaussian random number generator. Careful, h0 must be evaluated once at the initialisation of the ocean. Otherwise, at each time t, h will be evaluated from a different h~, and therefore will represent a different ocean.

We also have the following formula to compute the surface slopes.


To get the normal, we just compute : up_vector – slope

From there, we can code a heightmap generator to build our ocean.

Ocean heightmap. Grid of size 100x100. Computed in 77S.
Ocean heightmap. Grid of size 100×100. Computed in 77 seconds.

This gives an ocean surface quite realistic but it is not that satisfying yet. In real life, the ocean surface is not an heightmap

The ocean surface is not a heightmap !
The ocean surface is not a heightmap !

Therefore, Tessendorf presents the Choppy Waves algorithm. If it does not allow to create plunging waves as illustrated, it allow to creates sharp waves. This is very similar to Gerstner waves. Every point of the grid mesh is displaced on x and z to compress or dilate parts of the mesh; simply by adding D(x, t) to the x and z coordinates of the wave. This creates peak waves.


Heightmap ocean
Displaced grid ocean. Peak waves appear in the middle.


Leading to a final rendering

200×200 – 15 minutes to compute


Using a discrete Fourier Transform, we have the following code. However, as given in the caption of the pictures, computing the surface is very slow. This can be improved by using a Fast Fourier Transform, but this will be another story !

A small optimisation that can divide almost by 2 the computation times. For each h0~ that is computed, only part of them will impact strongly the final surface as some of them will lead to waves with amplitudes too small to be seen. What can be done in the initialisation therefore is to test the values of the h0~ and to keep only those above a thresold. Then, when computing h~, the sum will be done on much less k and will therefore run faster.

Grid of 100×100. No Thresholding (10 000 elements in the sum). Computed in 75 seconds


Grid of 100×100. Threshold at 0.01, leading to 5775 frequencies to compute the h~ sum. Computed in 44 seconds
class Ocean
    Ocean(float n, float m, 
          float lx, float lz, 
          Vector2 windir, float windspeed);
    void computeHeightmap(float time);
    Mesh* getMesh();
    int N, M; //Nombre de samples
    float Lx, Lz; //Distance entre deux samples
    Vector2 wind;
    Vector3* geometry;
    Vector3* normales;
    int nb_htilde;
    Complex* htilde0_tab;
    Complex* htilde0c_tab;
    Vector2* k_tab;
    unsigned int* indices;
    //h(x, t) + D(x, t) + n(x, t)
    void evaluate(Vector2 x, float t, 
                  float* h, Vector3* n, Vector2* disp);

    Complex htilde(Vector2 k, float t);    
    Complex htilde0(Vector2 k);
    double Phillips(Vector2 k);
    double dispersion(Vector2 k);
Ocean::Ocean(float n, float m, 
             float lx, float lz, 
             Vector2 windir, float windspeed)
    M = m;
    N = n;
    Lx = lx;
    Lz = lz;
    wind = windir.normalize() * windspeed;
    geometry = (Vector3*)malloc((M)*(N)*sizeof(Vector3));
    normales = (Vector3*)malloc((M)*(N)*sizeof(Vector3));
    htilde0_tab  = (Complex*)malloc((M)*(N)*sizeof(Complex));
    htilde0c_tab = (Complex*)malloc((M)*(N)*sizeof(Complex));
    k_tab = (Vector2*)malloc((M)*(N)*sizeof(Vector2));
    indices = (unsigned int*)malloc( 
                         (M-1)*(N-1)*sizeof(unsigned int)*6);

    nb_htilde = 0;    
    int index;

    for (int m = -M/2; m < M/2; m++)
    for (int n = -N/2; n < N/2; n++) {

        int n_prime = n + (N/2);
        int m_prime = m + (M/2);
        index = m_prime * (N) + n_prime;
        Vector2 k( 2 * M_PI * n / Lx,
        2 * M_PI * m / Lz);

        htilde0_tab[nb_htilde] = htilde0(k);
        if (htilde0_tab[nb_htilde].getModule() > 0.01)
            k_tab[nb_htilde] = k;
            htilde0c_tab[nb_htilde] = 

        geometry[index].x ( n * Lx / N );
        geometry[index].y ( 0.0f );
        geometry[index].z ( m * Lz / N );
        normales[index].x ( 0 );
        normales[index].y ( 1 );
        normales[index].z ( 0 );

    int indices_count = 0;
    for (int m_prime = 0; m_prime < M-1; m_prime++)
    for (int n_prime = 0; n_prime < N-1; n_prime++) {
        index = m_prime * (N) + n_prime;
        indices[indices_count++] = index; 
        indices[indices_count++] = index + 1;
        indices[indices_count++] = index + N;
        indices[indices_count++] = index + 1;
        indices[indices_count++] = index + N + 1;
        indices[indices_count++] = index + N;


/** To create and handle the meshes I use my own library that is not released **/
Mesh* Ocean::getMesh()
    Mesh* m = new Mesh();
    m->loadGeometry(geometry, (M)*(N));
    m->loadIndices(indices, (M-1)*(N-1)*6);
    m->addBuffer(M*N, 3, sizeof(Vector3), GL_FLOAT, normales);
    return m;

void Ocean::computeHeightmap(float time)
    int nb = 0;
    for(int n=-N/2; n<N/2; n+=1)
        //printf("%d / %d\n", nb, N);
        for(int m=-M/2; m<M/2; m+=1)
            int n_prime = n + (N/2);
            int m_prime = m + (M/2);
            int index = m_prime * (N) + n_prime;
            Vector2 x(n*Lx/N, m*Lz/M);
            float height;
            Vector3 normale;
            Vector2 disp;
            evaluate(x, time, &height, &normale, &disp);
            float lambda = 0.8;
                  ( n * Lx / N ) + lambda * disp.x() 
                  ( m * Lz / M ) + lambda * disp.y() 
            geometry[index].y( height );
            normales[index] = normale;
            normales[index].y( 1 );

void Ocean::evaluate(Vector2 x, float t, 
                     float* h, Vector3* norm_, Vector2* disp_)
    Complex height(0, 0);
    Vector3 normale;
    Vector2 disp;
    for(int i=0; i<nb_htilde; i++)
        Vector2 k = k_tab[i];
        if (k.length() < 0.000001)
        float kx = DotProduct(k, x);
        Complex h = htilde(k, t, i) * Complex(cos(kx), sin(kx));
        height += h;
        normale += Vector3(-k.x() * h.imaginary(), 
                           -k.y() * h.imaginary());
        disp += Vector2(-k.normalize().x() * h.imaginary(), 
                        -k.normalize().y() * h.imaginary());
    *h = height.real();
    *norm_ = normale;
    *disp_ = disp;

double Ocean::dispersion(Vector2 k)
    float w0 = 2.0f * M_PI / 200.0f;
    return floor(sqrt(9.8 * k.length()) / w0) * w0;

Complex Ocean::htilde(Vector2 k, float t, int idx)
    int n_prime = (k.x() * Lx) / (2* M_PI);
    n_prime += N/2;
    int m_prime = (k.y() * Lz) / (2* M_PI);
    m_prime += M/2;
    int index = idx;
    float omegat = dispersion(k) * t;
    float cos_ = cos(omegat);
    float sin_ = sin(omegat);
    Complex c0(cos_,  sin_);
    Complex c1(cos_, -sin_);
    return htilde0_tab[index] * c0 + htilde0c_tab[index] * c1;

Complex Ocean::htilde0(Vector2 k)
    double Ph = Phillips(k);
    Complex c = gaussianRandom();
    Complex r = c * sqrt(Ph/2.0);
    return r;

double Ocean::Phillips(Vector2 k)
    float V = wind.length();
    float L = (V*V) / 9.8;
    float A = 0.00001;

    float kl = k.length();
   if (kl < 0.000001)
        return 0.0;
    float kw = DotProduct(k.normalize(), wind.normalize());
    double ph = exp(-1 / (kl*L*kl*L) );
    ph /= kl*kl*kl*kl;
    ph *= kw*kw;
    ph *= A;
    float l = L * 0.001;
    ph *= exp(-kl*kl*l*l);
    return ph;

[1] A. Fournier and W. T. Reeves, “A simple model of ocean waves,” in ACM Siggraph Computer Graphics, 1986, vol. 20, pp. 75–84.
[2] G. A. Mastin, P. A. Watterberg, and J. F. Mareda, “Fourier synthesis of ocean scenes,” Computer Graphics and Applications, IEEE, vol. 7, no. 3, pp. 16–23, 1987.
[3] J. Tessendorf and others, “Simulating ocean water,” Simulating Nature: Realistic and Interactive Techniques. SIGGRAPH, vol. 1, 2001.

Fourier & co

Continuous Signal

What is the Fourier Transform ? This concept can be quite difficult to grasp at first. The Fourier Transform is a formula linking the spatial representation of a function with its representation in the frequency domain.

A signal s can be represented either by its value at a time t, s(t), or by its amplitude S for the frequency f, S(f). Example; we have a signal s

s(t) = A * cos( w* t )

Its Fourier Transform is equal to

S(f) =   √(π/2)  * [   A * δ (-w + f) + A * δ(w + f)  ]

Which pretty much means that it is worth 0 everywhere but when f = w. Which is normal. Since s is a cosine signal, it only expresses a single frequency. Therefore, when computing S(f) on every frequency, it only has an amplitude if f is the frequency of the cosine of s.

Usually, we consider that s and S are two formulation of the same function. To go back and forth between them, we use the following formulas


Keeping in mind that   b4e65f691d38a02a11091fac6285276b. You may note that those equations are linear. This means that is we set two functions s and s’ and that we note F(s) the Fourier transform of s, we have

F(s+s’) = F(s) + F(s’)
A * F(s) = F(A*s)

This leads a very nice property when convolving two functions; the Fourier transform of the convolution of a function s and a function s’ is actually equal to the product of their Fourier Transform.

F(s ∗ s’) = F(s)F(s’)

Discretely sampled data

The analytical functions stated above work well except in one case; when s is not defined analytically but rather as a set of samples.

Let’s do some sampling theory first. When sampling a data with a samplig rate fs , Nyquist theorem tells us that you con only recover frequency between fc and – f with

nyquist_limitThis has two consequences. First, if a function s(t) is bandwith limited to frequencies within the range [ -fc ,  fc ], then s is completely determined by its samples. The other  one is that if s(t) is NOT bandwith limited in [ -fc ,  fc ], then aliasing appear. Every frequency outside the range will be aliased (falsely translated). Once you have aliasing, there is pretty much nothing you can do.

To detect aliasing in a function, you can look at its fourier transform. If the function is correctly sampled, this transform should be 0 outside [ -fc ,  fc ]. Therefore, it should lean to 0 when the frequencies lean to fc . If it leans towards a finite value, there is great chances you have aliasing.

Once you have a correctly sampled signal s, how do you compute its fourier transform ? If you try to use the formulas given above, you will encounter a problem trying to integrate on values that are not given by the sampled. If we have N samples in input, we will be able to determine the Fourier coefficient for no more than N frequencies. We will then f estimate S(f) for


You may notice that the extreme values of n are the critical frequency range. From here, we replace the integrals with discrete sums


Fast Fourier Transform

Computing the discrete Fourier transform of N points appears to be a O(n²) problem. There is in fact a solution to compute it an a O(Nlog2(N)) process; this is the Fast Fourier Transform.

This algorithm starts with the statement that the Discrete Fourier transform of a set a samples N can be expresses ad the sum of the DFT of the even samples and the DFT of the odd samples. This is known as the Danielson-Laczos lemma.


And, interesting fact, it happens to be recursive ! However, this recursivity comes at the price that N is IMPERATIVELY a power of 2. If N is not a power of 2, you should pad it up with 0 until you get a power of 2.

We apply recursively this Lemma until we get to compute Fourier transforms of length one. This means that we have to compute a S from a single sample sn that has been successively even or odd log2(N) times. According to the equation of the discrete fourier transform, S is equal to its corresponding  sn.


Therefore, all we need to do is find which n correspond to the even/odd pattern. Fun fact, this is super easy to get ! Juste reverse the even/odd pattern and set e=0 and o=1. And there you get the binary value of n ! This is because all those even and odd classification are merely tests over the last bit of n.

The whole FFT algorithm starts by rearranging the N samples in bit reversed order; { 00, 01, 10, 11} becomes {00, 10, 01, 11} (you read the bit from right to left). Therefore, you can apply Danielson-Laczos lemma on the values pairs by pairs in this table until you get the final transform, after Nlog2(N) operations.


Ocean modelisation – Gerstner

Modeling ocean waves can be a difficult task. Waves are generated from the interaction between the wind and the ocean surface and this interaction is usually overly simplified during modelisation. There are two main classes of solution to generate waves, either using the Fourier spectra of the ocean [2], or representing the waves as a sum of trochoids [1].

This second representation is called the Gerstner model. When Fournier and Reeves presented this model, they wanted to be able to modelise waves breaking along the shores as well as those in the middle of the ocean. Here, we focus on deep ocean wave so we will put aside the changes induced by the shores.

Their model present several interesting caracteristics

  • Generate realistic waves
  • Generate realistic motion
  • Generate a wide range of waves, from ripples to high storm waves
  • Allow some user control

Although this model is realistic, it does NOT give any physical answer. This model might be useless to an oceanographer but it looks like the real ocean. In computer graphics and rendering, that is all we ask🙂

A small flaw with this paper would be that it only represent a single wavetrain. The rendered geometry is therefore lacking the real complexity of the ocean (which might be added by normal mapping though). However, it has the advantage that the whole ocean can be parametered only by a wind speed.

They use throchoïds to build their wave. This means that for each x point on the surface, this point is displaced on a circular orbit. This simulates the compression that happens in water motion, generating the waves.

Image take from [3]
In 1D with the y axis pointing up, the equations are

gerstner_1D_1Wwhere x is the position of the point and t the time. From this, any wavetrain can be generated depending on the parameters

  • h : the height of the wave
  • k : the spatial frequency of the wave
  • ω: the temporal frequency of the wave
Wave profile for h*k = 0.5 at t = 0
Wave profile for h*k = 1 at t = 0

But if you chose those parameters freely and randomly, you might see this kind of artefact appear

Wave profile for h * k = 1.5 at t = 0

This is the wave crossing itself. This kind of self intersection artefact appears when k*h > 1.

The dispersion relation of the waves gives us the relation between k and ω


In their paper, [1] propose several formula linking all the parameters to a wind speed V; allowing to generate a correct wavetrain with a single parameter.

  • h = 0.5 * 0.001 * 7.065 * V2.5
  • ω = (g/V) *  (2/3) 
  • k = ω² / g

where g is the gravitational constant.

Which gives as a primary result for deep ocean waves

V = 16
V = 12
V = 8
V = 4
float hfactor(float V)
    return 2 * 0.001 * 7.065 * pow(V, 2.5);

float wfactor(float V)
    return (9.8 / V) * sqrt(2.0/3.0);

float kfactor(float V)
    float w = wfactor(V);
    return (w*w)/9.8;

void main()
    vec3 dpos = position;
    float h = hfactor(windspeed);
    float k = kfactor(windspeed);
    float w = wfactor(windspeed);

    dpos.y = -h * cos(k*dpos.x - w*time);    
    dpos.x += h * sin(k*dpos.x - w*time);
    vertex_position = vec4(, 1);
    gl_Position = mvp * vertex_position;

This is nice but if we stick to [1], we are still limited to a single wavetrain. Also, in [1], the question of wind direction is handled as a rotation of the wavetrain but they do not get into the details.

Therefore, we will move to [4]. This paper from Tessendorf extends the model of [1] to handle several wavetrains. Furthermore, he handles the wind direction by simple dot products in the equations. The issue is that the relations between the variables presented above do not hold anymore as the equations in which they are to be applied changed. The only relation that hold is the dispersion relation.

The new equations are



  • x is the displaced point, and x0 its position at rest.
  • y is its coordinate in the up axis.
  • k is now a vector, called the wavevector. Its norm gives the wavelength and its direction gives the direction of the wave.
  • ω is the temporal frequency of the wave.
  • h is the amplitude of the wave
  • Φ is the phase of the wave

As you may note on those equations, what Tessendorf did was to sum several gerstner wavetrains with different parameters. The main issue I had while implementing this was to find a set of parameters that would give correct results.

In the end I went with generating random parameters. I used the following values :

  • ω = random number in [1:6]
  • |k| = ω² / 9.8
  • h = 1 / |k| * 0.05 (value chosen empirically).

almost_random_oceanBy using this and a fixed k, you can get a nice moving plane but nothing realistic. To add some randomness to it, I switched to the polar coordinates;

k (x, y) becomes k (|k|, θ). For x > 0 and y > 0, θ = arctan(y/x) ( ).

Then, I added a random angle from [-0.75, 0.75] to θ.

random_oceanHow does it translates in code;

What I did was to load a grid of points on the GPU and then, in a vertex shader, I displaced all the vertices according to the equations of Tessendorf. I generate the parameters for the waves in a pre process, and moves it to the GPU in a 1D Texture. The parameters to write are ω, h (name A in the code), |k|sin(θ) and |k|cos(θ).

CAREFUL, the texture must NOT be filtered🙂 ! You do not want the parameters from the first wavetrain to impact the second ! To avoid this use


and fill the texture with

unsigned int nbtrains = 128;
float* waves = (float*)malloc(sizeof(float)*4*nbtrains);
for(unsigned int i=0; i<nbtrains; i++)
        float w = (float)( rand()%30 + 7 ); // w in [7:36]
        w /= 7.0; //w in [1:36/7]
        float k_norm = (w*w)/9.8;
        float A = 1.0 / k_norm;
        A *= 0.05;//( (float)( rand()%100 ) / 100.0 );
        waves[i*4] = A;
        waves[i*4 + 1] = w;
        float theta_off = ( (float)(rand()%10 - 5)/7.0 );
        float x = winddir.x();
        float y = winddir.y();
        float theta = atan(y / x);
        theta += theta_off;
        float sin_theta = sin(theta);
        float cos_theta = cos(theta);
        waves[i*4 + 2] = k_norm * sin_theta;
        waves[i*4 + 3] = k_norm * cos_theta;

Then, all is left is to read this in the vertex shader.

 vec2 xz = dpos.xz;
 float h = 0;
 for(int i=0; i<nb_trains; i++)
        vec4 param = textureLod(waves_param, 
                   (float(i)/128.0) + 0.0001, 0);
        float A = param.x;
        float w = param.y;
        float sin_t = param.z;
        float cos_t = param.w;
        vec2 k = vec2(cos_t, sin_t);
        float phi = dot(k,dpos.xz) - w*time + i;
        xz -= normalize(k)*A*sin(phi);
        h += A * cos(phi);
 vertex_position = vec4(xz.x, h, xz.y, 1);

To get a nice ocean mesh🙂 !

Of course, with this solution the computations are looooooooooooooooooooooooooong …
There is a solution to increase the speed by using the FFT (Fast Fourier Transform) to compute the sum of sin/cos faster by I did not investigate it (yet)

[1] A. Fournier and W. T. Reeves, “A simple model of ocean waves,” in ACM Siggraph Computer Graphics, 1986, vol. 20, pp. 75–84.
[2] G. A. Mastin, P. A. Watterberg, and J. F. Mareda, “Fourier synthesis of ocean scenes,” Computer Graphics and Applications, IEEE, vol. 7, no. 3, pp. 16–23, 1987.
[3] E. Bruneton, F. Neyret, and N. Holzschuch, “Real-time Realistic Ocean Lighting using Seamless Transitions from Geometry to BRDF,” in Computer Graphics Forum, 2010, vol. 29, pp. 487–496.
[4] J. Tessendorf and others, “Simulating ocean water,” Simulating Nature: Realistic and Interactive Techniques. SIGGRAPH, vol. 1, 2001.

Filtering 101

Filtering ? Kesako ?

The process of rendering is based strongly on sampling, as an example, the scene to be rendered is sampled by the pixels of the final image. This implies that a correct rendering must respect the Nyquist limit; that the frequency of details in the scene must be 2 times lower than the sampling done by the pixels. Otherwise aliasing artefacts appear.

This aliasing can come from several reasons. It could be induced by the geometry if several triangles project onto the same pixel. This aliasing can also be induced by the spatial maps (textures, normal maps, horizon maps …) applied on the surface. To avoid this, GPU use Mipmapping, but this is not suitable in many cases.

The color of a pixel p is ideally the mean of all the visible contributions of the spatial maps in this pixel. This is DIFFERENT than the color of the visible mean contributions of the spatial maps.

The first teapot show the normal map at the original mipmapping level. The second shows the enlightenment obtained using just the mean normal obtained from mipmapping. The last teapot takes into account that several normals contributes to the final enlightenment of a point. Image taken from [1]
This becomes obvious with normal maps. Applying a normal map on an object, when seeing the object closely, the normal map is sampled correctly (first teapot). When seen from afar, if the normals are just mipmapped, the object becomes progressively smooth as solely the mean normal contributes to the enlightenment (second teapot). Instead, we would like to keep the information in each point that there is in fact a wide number of different normals contributing to the enlightenment, not only the mean one (third teapot).

This is the objective when filtering data. The point is the remove aliasing artefacts by finding new representations for the objects. We want to be able with a simple request (like a texel fetch or a single function evaluation) to retrieve all the information about the surface that projects into the pixel (in the case of normal maps, we want both a mean normal and the degree of roughness of the original surface).

Mathematical Background

This section is strongly inspired from Bruneton & Neyret [2]. I will try to detail the mathematical intuitions behind all this formalism.

Rendering is a integration problem. Integrating for each point over all the incoming light direction or integrating in a pixel over the whole patch of surface that projects onto it. Here, we will focus on this second integration.

Here, we render an object S. This object is represented in memory with a coarse geometry (A) and several spatial maps. We want to compute the value I of the pixel onto which A projects. This is the integral over S of the local illumination equation (in our case, this equations returns the color of a pixel depending on the light direction, camera direction and the spatial maps). Picture from [2].
To solve an integral on a computer, several methods exists that does numerical integration. They usually sum all the samples of the function over the integration domain. This is quite slow and returns just an approximation of the integral value. Therefore, we would rather to be able to solve our integral analytically; expressing it with a function g that can be evaluated on the computer.

In our case, we rely on spatial maps. Thus, our function will be parametrized by those maps at different level of resolutions, depending on S. Once we found our function g, we will use textures to store its parameters and mipmapping to get them at the right resolution. This implies that all those parameters must also combine linearly (that the averaging done by the mipmapping process gives a correct result).

For example, if we put a function f parametrized by a map M. The color of the pixel I is computed with

We want to re express the values in M to get a new map M’  such as



This leaves the question of the function g. Most of the time, g is based on the probability distribution of the M values over S. Often, it is a Gaussian function, parametrized only by its mean and its variance. The mean value combines linearly and can be stored in a texture. The variance does not combine linearly. Instead, we use the moments of the function and exploit the fact that the variance can be compute from the mean value and the mean square value. It is then stored in a second texture, containing the squared values of the function.


It might be easier to understand all this with a concrete practical case

We have a texture with gaussianly distributed values (we can check this by looking at the histogram of the image to see if it looks like a gaussian )


Applying this texture on an infinite plane, we want to color every point above a value t with a color Ch and all other points with a color Cl.

The rendering equation therefore becomes


Where H is the heaviside function (H(x) = 0 if x < 0, 1 otherwise), h(x) is the height at a point x and S is the surface patch. We keep going on with the math, seperating the integrals from Ch and Cl that are constants.step1

Now, we change the integrals. Instead of integrating over S, we integrate over D, the repartition domain of the heights. Instead of summing up all the heights into the pixel, we sum up all the possible heights, weighted by the probability they appear in the pixel.


NDF(S, h) is the probability that h is present on the patch S.

Then, by using the fact that H(x) = 0 when x < 0, we reduce the integration domain. As H(x) = 1 when x > 0, we can also remove it from the formula (we suppose that D = [0, 1])


Which is equal to


Where CDF_h(S, t) is the cumulative distribution function of h over the patch S until the value t (the CDF of a probability distribution function is by definition, an analytical function of its integrand).

All those maths means that to compute the filtered color of our pixel, we need to be able to evaluate for every S the CDF of the distribution function. Here, our heights are gaussianly distributed. Thus, the CDF is the CDF of a gaussian function


So, to filter our pixel, we only need the be able to get the mu and sigma in each pixel. And this can be done with textures and mipmapping.

Loading the noise tile as a texture and mipmapping already gives us mu. The only trick is to get the sigma value. To get it, we use this property

variance_esperancewhere Var(x) = sigma and E[X] is the expected value of x (here, it is its mean, and therefore it is the value computed by the mipmapping process). E[X] is the mu value given by the texture containing the tile, and E[X²] can be obtained by sampling a second texture containing the squared noise values.

In code terms it means

float t = 0.47;

float height = texture(heightmap, tex_coord).r;
float height2 = texture(heightmap_squared, tex_coord).r;
float mu = height;
float sigma2 = max(1e-4, height2 - height*height);
float max_cdf = 0.5 * ( 1 +  erf((1 -mu)/sqrt(2.0 * sigma2)));
float cdf = 0.5 * (1 + erf((t -mu)/sqrt(2.0 * sigma2)));
color = vec4( (max_cdf - cdf) * vec3(0.7, 0.7, 0.7) + cdf * vec3(0, 0.5, 0), 1);

And in picture terms

Without filtering
With filtering


[1] M. Olano and D. Baker, “LEAN mapping,” in Proceedings of the 2010 ACM SIGGRAPH symposium on Interactive 3D Graphics and Games, 2010, pp. 181–188.

[2] E. Bruneton and F. Neyret, “A survey of nonlinear prefiltering methods for efficient and accurate surface shading,” Visualization and Computer Graphics, IEEE Transactions on, vol. 18, no. 2, pp. 242–260, 2012.