Blog of roxlu, co-founder of Apollo Media. Contact info[shift+2]apollomedia.nl.

Height Field Simulation on GPU

The code below shows a first rough implementation of a height field simulation on GPU. It was an experiment to test different approaches to perform calculations on the GPU; this one uses plain textures from which we read the current values and another set of textures into which the new results are written.

The height field algorithm diffuses a certain value In my case I used this value as height. On the CPU this looks something like the code below. We need two arrays, one array float u0[FIELD_NN] that contains the current heights, and another one float u1[FIELD_NN] into which we write the changed values using the current velocity array (float v[FIELD_NN]). After traversing all the cells in the grid, we copy the new values from u1 into u0 so the newest values will be used during the next iteration.

The value of c needs to be in a certain range if you want to get good results, see the height field slides from the link above.

Download the code here. Tested only on a mac 10.9.

Height Field on CPU

void Water::updateHeightField(float dt) {
  float c = 1.0;
  float max_c = (1.0 / dt);
  if(c > max_c) {
    printf("Warning: invalid C value\n");
    return;
  }
 
  dt = dt * 4.9;
  for(int j = 1; j < FIELD_N - 1; ++j) {
    for(int i = 1; i < FIELD_N - 1; ++i) {
      int current = ((j + 0) * FIELD_N) + (i + 0);
      int right   = ((j + 0) * FIELD_N) + (i + 1);
      int left    = ((j + 0) * FIELD_N) + (i - 1);
      int top     = ((j + 1) * FIELD_N) + (i + 0);
      int bottom  = ((j - 1) * FIELD_N) + (i + 0);
      float f = ((c*c) * ((u0[right] + u0[left] + u0[bottom] + u0[top])) - (4.0 * u0[current])) / 4.0 ;
      if(f > 0.1) {
        f = 0.1;
      }
      else if(f < -0.1) {
        f = -0.1;
      }
 
      v[current] = v[current] + f * dt;
      u1[current] = u0[current] + v[current] * dt;
      v[current] *= 0.995;
 
    }
  }
  memcpy(u0, u1, sizeof(u0));
}

On the GPU, we need to ping/pong between the source values and the destination values for both the height values and the velocity values, so we need 4 textures for those two values. In the code below I also have a texture into which we write the normals that we calculate in a separate step.

HeightField.h

/*
 
  HeightField
  -----------
 
  This class implements a basic height field diffusion as described in 
  http://www.matthiasmueller.info/talks/gdc2008.pdf
 
  We perform the following steps:
 
  - Diffuse the current height values according to some
    velocity. We use u0 and u1 as the height values. We need two
    because we're ping/ponging the values. The current velocity
    at which the height field diffuses is stored in v0 and v1.
    After each time step we change `state_diffuse` which toggles
    the reading/writing from u0/v0 and u1/v1.
 
  - Once we've done the diffuse step, we perform a processing step
    where we calculate the world position of the heightfield. These
    height values are stored in `tex_out_pos`.
 
  - When we have the positions in `tex_out_pos` we perform on more
    step where we calculate the normals, and values we might need for
    water rendering.
 
  - To render the height field you can use `vertices_vbo` and `vertices_vao`.
    They are setup in such a way that you'll get one in attribute in your shader
    called `tex` which contains the row/column into the position texture. Use this
    to set gl_Position.
 
  Things to know / improve:
 
  - We're not reusing the attribute less GL_VERTEX_SHADER;
 
  - You can use the vertices_vao/vbo to render the height field, see the 
    debug shaders.
 
 */
#ifndef ROXLU_HEIGHT_FIELD_H
#define ROXLU_HEIGHT_FIELD_H
 
#define ROXLU_USE_ALL
#include <tinylib.h>
#include <vector>
 
struct HeightFieldVertex {
  HeightFieldVertex(vec2 texcoord):tex(texcoord){}
  HeightFieldVertex(){}
  vec2 tex;    
};
 
class HeightField {
 
 public:
  HeightField();
  bool setup(int w, int h);
  void update(float dt);     /* diffuses the height field */
  void process();            /* processes the current values, calculates normals, create position texture etc.. */
  void debugDraw();
 
  void print();               /* print some debug info */
 
 private: 
  bool setupDiffusing();     /* setup GL state for the diffusion step */
  bool setupProcessing();    /* setup GL state for the processing step; calculates normals, positions, texcoord etc.. using the current field values */
  bool setupDebug();         /* setup GL state for debugging */
  bool setupVertices();      /* create the triangle mesh (or the order of vertices, position is extracted from the position texture) */
 
 public:
 
  /* diffusion of height field */
  GLuint field_fbo;
  GLuint tex_u0;              /* height value */
  GLuint tex_u1;              /* height value */
  GLuint tex_v0;              /* velocity at which u diffuses */
  GLuint tex_v1;              /* velocity at which u diffuses */
  GLint u_dt;                 /* reference to our dt uniform */
  GLuint field_vao;           /* we need a vao to render attribute less vertices */
  Program field_prog;         /* this program does the diffuse step */
  int state_diffuse;          /* toggles between 0 and 1 to ping/pong the diffuse/vel textures */
 
  /* general info */
  int field_size;             /* the size of our rectangular height field */
  int win_w;                  /* window width, use to reset the viewport */
  int win_h;                  /* window height, used to reset the viewport */
 
  /* used to process the height field and extract some usefull data */
  GLuint process_fbo;         /* we use a separate FBO to perform the processing step so we have some space for extra attachments */
  GLuint tex_out_norm;        /* the GL_RGB32F texture that will keep our calculated normals */
  GLuint tex_out_pos;         /* the GL_RGB32F texture that will keep our positions */
  GLuint tex_out_texcoord;    /* the GL_RG32F texture that will keep our texcoords */
  Program process_prog;       /* the program we use to calculate things like normals, etc.. */
  Program pos_prog;           /* the program we use to calculate the positions */
 
  /* used to debug draw */
  Program debug_prog;         /* debug shaders, shows how to set gl_Position */
  mat4 pm;                    /* projection matrix */
  mat4 vm;                    /* view matrix */
 
  /* vertices */
  std::vector<HeightFieldVertex> vertices;    /* The vertices that you can use to render a triangular height field, see the debug shader */
  GLuint vertices_vbo;                        /* VBO that holds the HeightFieldVertex data that forms the height field triangular grid */
  GLuint vertices_vao;                        /* The VAO to draw the height field vertices */
};
 
#endif

HeightField.cpp

#include "HeightField.h"
 
HeightField::HeightField() 
  :field_fbo(0)
  ,tex_u0(0)
  ,tex_u1(0)
  ,tex_v0(0)
  ,tex_v1(0)
  ,field_size(128)
  ,field_vao(0)
  ,state_diffuse(0)
  ,win_w(0)
  ,win_h(0)
  ,process_fbo(0)
  ,tex_out_norm(0)
  ,tex_out_pos(0)
  ,tex_out_texcoord(0)
  ,vertices_vbo(0)
  ,vertices_vao(0)
{
}
 
bool HeightField::setup(int w, int h) {
 
  if(!w || !h) {
    printf("Error: invalid width/height: %d x %d\n", w, h);
    return false;
  }
 
  win_w = w;
  win_h = h;
 
  glGenVertexArrays(1, &field_vao);
 
  if(!setupDiffusing()) {
    printf("Error: cannot set GL state for the diffuse step.\n");
    return false;
  }
 
  if(!setupProcessing()) {
    printf("Error: cannot setup the GL state for the processing.\n");
    return false;
  }
 
  if(!setupVertices()) {
    printf("Error: cannot setup the vertices for the height field.\n");
    return false;
  }
 
  if(!setupDebug()) {
    printf("Error: cannot setup the GL state for debugging.\n");
    return false;
  }
 
  return true;
}
 
bool HeightField::setupVertices() {
  glGenVertexArrays(1, &vertices_vao);
  glBindVertexArray(vertices_vao);
 
  glGenBuffers(1, &vertices_vbo);
  glBindBuffer(GL_ARRAY_BUFFER, vertices_vbo);
 
  std::vector<HeightFieldVertex> tmp(field_size * field_size, HeightFieldVertex());
  for(int j = 0; j < field_size; ++j) {
    for(int i = 0; i < field_size; ++i) {
      int dx = j * field_size + i;
      tmp[dx].tex.set(i, j);
    }
  }
 
  for(int j = 0; j < field_size-1; ++j) {
    for(int i = 0; i < field_size-1; ++i) {
      int a = (j + 0) * field_size + (i + 0);
      int b = (j + 0) * field_size + (i + 1);
      int c = (j + 1) * field_size + (i + 1);
      int d = (j + 1) * field_size + (i + 0);
      vertices.push_back(tmp[a]);
      vertices.push_back(tmp[b]);
      vertices.push_back(tmp[c]);
 
      vertices.push_back(tmp[a]);
      vertices.push_back(tmp[c]);
      vertices.push_back(tmp[d]);
    }
  }
 
  glBufferData(GL_ARRAY_BUFFER, sizeof(HeightFieldVertex) * vertices.size(), vertices[0].tex.ptr(), GL_STATIC_DRAW);
 
  glEnableVertexAttribArray(0); // tex
  glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, sizeof(HeightFieldVertex), (GLvoid*)0);
  return true;
}
 
bool HeightField::setupDebug() {
 
  pm.perspective(60.0f, float(win_w)/win_h, 0.01f, 100.0f);
  //vm.lookAt(vec3(0.0f, 35.0f, 0.0f), vec3(0.0f, 0.0f, 0.1f), vec3(0.0f, 1.0f, 0.0f));
  vm.translate(0.0f, 0.0f, -30.0f);
  vm.rotateX(30 * DEG_TO_RAD);
 
  const char* atts[] = { "a_tex" };
  debug_prog.create(GL_VERTEX_SHADER, rx_to_data_path("height_field_debug.vert"));
  debug_prog.create(GL_FRAGMENT_SHADER, rx_to_data_path("height_field_debug.frag"));
  debug_prog.link(1, atts);
 
  glUseProgram(debug_prog.id);
  glUniformMatrix4fv(glGetUniformLocation(debug_prog.id, "u_pm"), 1, GL_FALSE, pm.ptr());
  glUniformMatrix4fv(glGetUniformLocation(debug_prog.id, "u_vm"), 1, GL_FALSE, vm.ptr());
  glUniform1i(glGetUniformLocation(debug_prog.id, "u_pos_tex"), 0);
 
  return true;
}
 
bool HeightField::setupProcessing() {
 
  glGenTextures(1, &tex_out_norm);
  glBindTexture(GL_TEXTURE_2D, tex_out_norm);
  glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB32F, field_size, field_size, 0, GL_RGB, GL_FLOAT, 0);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
 
  glGenTextures(1, &tex_out_pos);
  glBindTexture(GL_TEXTURE_2D, tex_out_pos);
  glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB32F,field_size, field_size, 0, GL_RGB, GL_FLOAT, 0);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
 
  glGenTextures(1, &tex_out_texcoord);
  glBindTexture(GL_TEXTURE_2D, tex_out_texcoord);
  glTexImage2D(GL_TEXTURE_2D, 0, GL_RG32F, field_size, field_size, 0, GL_RG, GL_FLOAT, 0);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
 
  glGenFramebuffers(1, &process_fbo);
  glBindFramebuffer(GL_FRAMEBUFFER, process_fbo);
 
  glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_2D, tex_out_pos, 0);
  glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT1, GL_TEXTURE_2D, tex_out_norm, 0);
  glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT2, GL_TEXTURE_2D, tex_out_texcoord, 0);
 
  if(glCheckFramebufferStatus(GL_FRAMEBUFFER) != GL_FRAMEBUFFER_COMPLETE) {
    printf("Error: process framebuffer not complete.\n");
    return false;
  }
 
  // Position processing
  const char* pos_frags[] = { "out_pos" };
  pos_prog.create(GL_VERTEX_SHADER, rx_to_data_path("height_field.vert"));
  pos_prog.create(GL_FRAGMENT_SHADER, rx_to_data_path("height_field_pos.frag"));
  pos_prog.link(0, NULL, 1, pos_frags);
  glUseProgram(pos_prog.id);
  glUniform1i(glGetUniformLocation(pos_prog.id, "u_height_tex"), 0);
  glUniform1i(glGetUniformLocation(pos_prog.id, "u_vel_tex"), 1);
 
 
  // Extra processing
  const char* process_frags[] = { "out_norm", "out_tex" };
  process_prog.create(GL_VERTEX_SHADER, rx_to_data_path("height_field.vert"));
  process_prog.create(GL_FRAGMENT_SHADER, rx_to_data_path("height_field_process.frag"));
  process_prog.link(0, NULL, 2, process_frags);
  glUseProgram(process_prog.id);
  glUniform1i(glGetUniformLocation(process_prog.id, "u_height_tex"), 0);
  glUniform1i(glGetUniformLocation(process_prog.id, "u_vel_tex"), 1);
  glUniform1i(glGetUniformLocation(process_prog.id, "u_pos_tex"), 2);
 
  glBindFramebuffer(GL_FRAMEBUFFER, 0);
 
  return true;
}
 
bool HeightField::setupDiffusing() {
 
  // some text data
  float* u = new float[field_size * field_size];
  float* v = new float[field_size * field_size];
  int upper = 50;
  int lower = 30;
  for(int j = 0; j < field_size; ++j) {
    for(int i = 0; i < field_size; ++i) {
      u[j * field_size + i] = 0.0f;
      v[j * field_size + i] = 0.0f;
      if(i > lower && i < upper && j > lower && j < upper) {
        u[j * field_size + i] = 3.5;
      }
    }
  }
 
  glGenTextures(1, &tex_u0);
  glBindTexture(GL_TEXTURE_2D, tex_u0);
  glTexImage2D(GL_TEXTURE_2D, 0, GL_R32F, field_size, field_size, 0, GL_RED, GL_FLOAT, u);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
 
  glGenTextures(1, &tex_u1);
  glBindTexture(GL_TEXTURE_2D, tex_u1);
  glTexImage2D(GL_TEXTURE_2D, 0, GL_R32F, field_size, field_size, 0, GL_RED, GL_FLOAT, u);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
 
  glGenTextures(1, &tex_v0);
  glBindTexture(GL_TEXTURE_2D, tex_v0);
  glTexImage2D(GL_TEXTURE_2D, 0, GL_R32F, field_size, field_size, 0, GL_RED, GL_FLOAT, v);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
 
  glGenTextures(1, &tex_v1);
  glBindTexture(GL_TEXTURE_2D, tex_v1);
  glTexImage2D(GL_TEXTURE_2D, 0, GL_R32F, field_size, field_size, 0, GL_RED, GL_FLOAT, v);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
 
 
  glGenFramebuffers(1, &field_fbo);
  glBindFramebuffer(GL_FRAMEBUFFER, field_fbo);
  glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_2D, tex_u0, 0);
  glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT1, GL_TEXTURE_2D, tex_u1, 0);
  glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT2, GL_TEXTURE_2D, tex_v0, 0);
  glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT3, GL_TEXTURE_2D, tex_v1, 0);
 
  if(glCheckFramebufferStatus(GL_FRAMEBUFFER) != GL_FRAMEBUFFER_COMPLETE) {
    printf("Error: diffuse framebuffer not complete.\n");
    return false;
  }
 
  const char* frags[] = { "out_height", "out_vel" } ;
  field_prog.create(GL_VERTEX_SHADER, rx_to_data_path("height_field.vert"));
  field_prog.create(GL_FRAGMENT_SHADER, rx_to_data_path("height_field.frag"));
  field_prog.link(0, NULL, 2, frags);
 
  glUseProgram(field_prog.id);
  glUniform1i(glGetUniformLocation(field_prog.id, "u_height_tex"), 0);
  glUniform1i(glGetUniformLocation(field_prog.id, "u_vel_tex"), 1);
  u_dt = glGetUniformLocation(field_prog.id, "u_dt");
 
  glBindFramebuffer(GL_FRAMEBUFFER, 0);
 
  delete[] v;
  delete[] u; 
  u = NULL;
  v = NULL;
  return true;
}
 
void HeightField::update(float dt) {
  glDisable(GL_DEPTH_TEST);
  glDisable(GL_BLEND);
 
  glViewport(0, 0, field_size, field_size);
  glBindFramebuffer(GL_FRAMEBUFFER, field_fbo);
  glUseProgram(field_prog.id);
  glBindVertexArray(field_vao);
  glUniform1f(u_dt, dt);
 
  state_diffuse = 1 - state_diffuse;
 
  if(state_diffuse == 0) {
    // read from u0, write to u1
    // read from v0, write to v1
    GLenum drawbufs[] = { GL_COLOR_ATTACHMENT1, GL_COLOR_ATTACHMENT3 };
    glDrawBuffers(1, drawbufs);
 
    glActiveTexture(GL_TEXTURE0);  glBindTexture(GL_TEXTURE_2D, tex_u0);
    glActiveTexture(GL_TEXTURE2);  glBindTexture(GL_TEXTURE_2D, tex_v0);
 
    glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
  }
  else {
    // read from u1, write to u0
    // read from v1, write to v0
    GLenum drawbufs[] = { GL_COLOR_ATTACHMENT0, GL_COLOR_ATTACHMENT2 };
    glDrawBuffers(2, drawbufs);
 
    glActiveTexture(GL_TEXTURE0);  glBindTexture(GL_TEXTURE_2D, tex_u1);
    glActiveTexture(GL_TEXTURE3);  glBindTexture(GL_TEXTURE_2D, tex_v1);
    glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
  }
 
  glBindFramebuffer(GL_FRAMEBUFFER, 0);
  glViewport(0, 0, win_w, win_h);
}
 
void HeightField::process() {
  glBindFramebuffer(GL_FRAMEBUFFER, process_fbo);
  glViewport(0, 0, field_size, field_size);
  glBindVertexArray(field_vao);
 
  glActiveTexture(GL_TEXTURE0); glBindTexture(GL_TEXTURE_2D, tex_u0);
  glActiveTexture(GL_TEXTURE1); glBindTexture(GL_TEXTURE_2D, tex_v0);
 
  {
    // Calculate positions.
    glUseProgram(pos_prog.id);
 
    GLenum drawbufs[] = { GL_COLOR_ATTACHMENT0 } ; // , GL_COLOR_ATTACHMENT1, GL_COLOR_ATTACHMENT2 } ;
    glDrawBuffers(1, drawbufs);
    glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
  }
 
  {
    // Use positions to calc normals, etc..
    glUseProgram(process_prog.id);
    GLenum drawbufs[] = { GL_COLOR_ATTACHMENT1, GL_COLOR_ATTACHMENT2 };
    glDrawBuffers(2, drawbufs);
 
    glActiveTexture(GL_TEXTURE2); 
    glBindTexture(GL_TEXTURE_2D, tex_out_pos);
 
    glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
  }
 
  glBindFramebuffer(GL_FRAMEBUFFER, 0);
}
 
void HeightField::debugDraw() {
  glViewport(0, 0, win_w, win_h);
  glBindFramebuffer(GL_FRAMEBUFFER, 0); // tmp
  glDisable(GL_DEPTH_TEST);
 
  glBindVertexArray(vertices_vao);
  glUseProgram(debug_prog.id);
 
  glActiveTexture(GL_TEXTURE0);
  glBindTexture(GL_TEXTURE_2D, tex_out_pos);
 
  glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
  //glDrawArrays(GL_TRIANGLES, 0, vertices.size());
  glDrawArrays(GL_POINTS, 0, vertices.size());
  glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
 
#if 1
  glBindFramebuffer(GL_READ_FRAMEBUFFER, field_fbo);
  glBindFramebuffer(GL_DRAW_FRAMEBUFFER, 0);
  glReadBuffer(GL_COLOR_ATTACHMENT0);
  glBlitFramebuffer(0, 0, field_size, field_size, 0, 0, field_size, field_size, GL_COLOR_BUFFER_BIT, GL_LINEAR);
#endif
}
 
 
void HeightField::print() {
  printf("heightfield.tex_u0: %d\n", tex_u0);
  printf("heightfield.tex_u1: %d\n", tex_u1);
  printf("heightfield.tex_v0: %d\n", tex_v0);
  printf("heightfield.tex_v1: %d\n", tex_v1);
  printf("heightfield.tex_out_norm: %d\n", tex_out_norm);
  printf("heightfield.tex_out_texcoord: %d\n", tex_out_texcoord);
  printf("heightfield.tex_out_pos: %d\n", tex_out_pos);
}

height_field.vert

#version 150
 
const vec2[] pos = vec2[4](
  vec2(-1.0, 1.0),
  vec2(-1.0, -1.0),
  vec2(1.0, 1.0),
  vec2(1.0, -1.0)
);
 
const vec2[] tex = vec2[4](
  vec2(0.0, 1.0),
  vec2(0.0, 0.0),
  vec2(1.0, 1.0),
  vec2(1.0, 0.0)
);
 
out vec2 v_tex;
 
void main() {
  gl_Position = vec4(pos[gl_VertexID], 0.0, 1.0);
  v_tex = tex[gl_VertexID];
}

height_field.frag

/*
 
  Height Field Water Simulation
  -----------------------------
  - See http://www.matthiasmueller.info/talks/gdc2008.pdf
  - The value of c must be limited to: 1.0/dt, so 30fps, 1.0/0.033 = 33.33.
 
 */
 
#version 150
 
uniform sampler2D u_height_tex;
uniform sampler2D u_vel_tex;
uniform float u_dt;
out float out_height;
out float out_vel;
in vec2 v_tex;
 
void main() {
 
  float u        = texture(u_height_tex, v_tex).r;
  float u_right  = textureOffset(u_height_tex, v_tex, ivec2( 1.0,  0.0)).r;
  float u_left   = textureOffset(u_height_tex, v_tex, ivec2(-1.0,  0.0)).r;
  float u_top    = textureOffset(u_height_tex, v_tex, ivec2( 0.0, -1.0)).r;
  float u_bottom = textureOffset(u_height_tex, v_tex, ivec2( 0.0,  1.0)).r;
 
  float c = 50.0;
  float f = ((c*c) * ((u_right + u_left + u_top + u_bottom) - (4.0 * u)) ) / 4.0;
 
  float v = texture(u_vel_tex, v_tex).r + f * u_dt;
  out_height = u + v * u_dt;
  out_vel = v * 0.996;
}

height_field_debug.vert

#version 150
 
in vec2 a_tex;
uniform mat4 u_pm;
uniform mat4 u_vm;
uniform sampler2D u_pos_tex;
 
void main() {
  vec4 world_pos = vec4(texelFetch(u_pos_tex, ivec2(a_tex), 0).rgb, 1.0);
  gl_Position = u_pm * u_vm * world_pos;
}

height_field_debug.frag

#version 150
 
out vec4 fragcolor;
 
void main() {
  fragcolor = vec4(1.0, 0.0, 0.0, 1.0);
}

height_field_pos.frag

#version 150
 
uniform sampler2D u_height_tex;
uniform sampler2D u_vel_tex;
out vec3 out_pos;
in vec2 v_tex;
 
const float size_x = 35.0;
const float size_y = 35.0;
const float field_size = 128.0;
const float half_x = size_x * 0.5;
const float half_y = size_y * 0.5;
 
void main() {
  float u = texture(u_height_tex,v_tex).r;
  out_pos = vec3(-half_x + (v_tex.s * size_x), 
                  u, 
                 -half_y + (v_tex.t * size_y));
}

height_field_process.frag

/*
 
  - Uses the calculated position to calculate the height values.
  - I'm not sure is I need to calculate dx/dy using a stepsize of 2 or 1
 
 */
#version 150
 
uniform sampler2D u_pos_tex;
uniform sampler2D u_height_tex;
uniform sampler2D u_vel_tex;
out vec3 out_norm;
out vec2 out_tex;
in vec2 v_tex;
 
void main() {
  vec3 center = texture(u_pos_tex, v_tex).rgb;
  vec3 right  = textureOffset(u_pos_tex, v_tex, ivec2(1,  0)).rgb;
  vec3 top    = textureOffset(u_pos_tex, v_tex, ivec2(0, -1)).rgb;
 
#if 1
  vec3 dx = right - center;
  vec3 dy = top - center;
#else
  vec3 dx = center - right;
  vec3 dy = center - top;
#endif
 
  out_norm = normalize(cross(dx, dy));
  out_tex = v_tex;
}
 
//vec3 left   = textureOffset(u_pos_tex, v_tex, vec2(-1.0,  0.0)).rgb;
//vec3 bottom = textureOffset(u_pos_tex, v_tex, vec2( 0.0,  1.0)).rgb;
//vec3 dx = left - right;
//vec3 dy = top - bottom;