clouds.cpp

   1 #include "debug.hpp"
   2 #include "program.hpp"
   3 #include "simulation.h"
   4 #include <GL/glew.h>
   5 #include <GLUT/glut.h>
   6 #include <array>
   7 #include <cmath>
   8 #include <cstdio>
   9 #include <cstdlib>
  10 #include <chrono>
  11 #include <glm/ext.hpp>
  12 #include <glm/glm.hpp>
  13 #include <sys/stat.h>
  14 #include <vector>
  15
  16 #pragma clang diagnostic ignored "-Wdeprecated-declarations"
  17
  18 enum Mode { render, debugContDist, debugColor, debugProbExt, debugProbAct };
  19 Mode curMode = render;
  20
  21 using namespace std;
  22 using namespace glm;
  23
  24 const float metaballR = 1.f / 16.f;
  25 inline float metaballField(float r) {
  26   if (r > 1)
  27     return 0;
  28   const float a = r / (1);
  29   return (-4.f / 9.f * powf(a, 6)) + (17.f / 9.f * powf(a, 4)) -
  30          (22.f / 9.f * powf(a, 2)) + 1;
  31 }
  32
  33 const float normalizationFactor = (748.f / 405.f) * M_PI;
  34
  35 void checkError() {
  36   if (GLenum e = glGetError()) {
  37     fprintf(stderr, "%s\n", gluErrorString(e));
  38     abort();
  39   }
  40 }
  41
  42 GLuint bbProg, sunProg;
  43 GLuint bbVao;
  44
  45 // Here we need to generate n_q textures for different densities of metaballs
  46 // These textures then go on the billboards
  47 // The texture stores attenuation ratio?
  48
  49 #define NQ 64
  50 GLuint bbTexIds[NQ];
  51
  52 // Stores attenuation ratio inside r channel
  53 // Should be highest value at center
  54 void precalculateBillboardTextures() {
  55   fprintf(stderr, "Calculating billboard textures...\n");
  56   glGenTextures(NQ, bbTexIds);
  57
  58   for (int d = 0; d < NQ; d++) {
  59     float data[32 * 32];
  60     for (int j = 0; j < 32; j++) {
  61       for (int i = 0; i < 32; i++) {
  62         // TODO: properly calculate this instead of whatever this is
  63         float r = distance(vec2(i, j), vec2(16, 16)) / 16;
  64         float density = (float)d / NQ;
  65         data[i + j * 32] =
  66             1 - (3 * density * 0.7 * (metaballField(r) / normalizationFactor));
  67       }
  68     }
  69
  70     mkdir("bbtex", 0777);
  71     char path[32];
  72     snprintf(path, 32, "bbtex/%i.tga", d);
  73     saveGrayscale(data, 32, 32, path);
  74
  75     glBindTexture(GL_TEXTURE_2D, bbTexIds[d]);
  76
  77     glTexImage2D(GL_TEXTURE_2D, 0, GL_RED, 32, 32, 0, GL_RED, GL_FLOAT, data);
  78     glGenerateMipmap(GL_TEXTURE_2D); // required, otherwise texture is blank
  79
  80     fprintf(stderr, "\r%i out of %i densities calculated%s", d + 1, NQ,
  81             d == NQ - 1 ? "\n" : "");
  82   }
  83 }
  84
  85 struct Metaball {
  86   vec3 pos;
  87   ivec3 coords;
  88   /** Density */
  89   float d;
  90   vec4 col;
  91 };
  92
  93 array<Metaball, CLOUD_DIM_X * CLOUD_DIM_Y * CLOUD_DIM_Z> metaballs;
  94
  95 const float cloudScale = metaballR;
  96 const float metaballScale = metaballR * 3.f;
  97 Clouds cs;
  98
  99 void calculateMetaballs() {
 100   stepClouds(&cs);
 101   /* for (int i = 0; i < 256; i++) { */
 102   /*   float x = ((float)rand()/(float)(RAND_MAX) - 0.5) * 2; */
 103   /*   float y = ((float)rand()/(float)(RAND_MAX) - 0.5) * 2; */
 104   /*   float z = ((float)rand()/(float)(RAND_MAX) - 0.5) * 2; */
 105   /*   float r = (float)rand()/(float)(RAND_MAX) * 1; */
 106   /*   Metaball m = {{x,y,z}, r}; */
 107   /*   metaballs.push_back(m); */
 108   /* } */
 109   for (int i = 0; i < CLOUD_DIM_X; i++) {
 110     for (int j = 0; j < CLOUD_DIM_Y; j++) {
 111       for (int k = 0; k < CLOUD_DIM_Z; k++) {
 112         Metaball m = {vec3(i, j, k) * vec3(cloudScale), {i, j, k}};
 113         /* m.pos = (m.pos * vec3(2)) - (cloudScale / 2); */
 114         m.pos -= vec3(CLOUD_DIM_X, CLOUD_DIM_Y, CLOUD_DIM_Z) * cloudScale / 2.f;
 115         m.d = cs.q[i][j][k];
 116         /* m.d = 0; */
 117         metaballs[i * CLOUD_DIM_Y * CLOUD_DIM_Z + j * CLOUD_DIM_Z + k] = m;
 118       }
 119     }
 120   }
 121   /* for (int z = 0; z < CLOUD_DIM_Z; z++) */
 122   /*   metaballs[32 * CLOUD_DIM_Y * CLOUD_DIM_Z + 32 * CLOUD_DIM_Z + z].d = 1;
 123    */
 124 }
 125
 126 vec3 sunPos = {0, 10, 0}, sunDir = {0, -1, 0};
 127 size_t envColorIdx = 0;
 128 // First color is sun color, second is sky color
 129 std::array<std::array<vec4, 2>, 3> envColors{
 130     {{vec4(1, 1, 1, 1), vec4(0.9, 1, 1, 1)},
 131      {vec4(0.939, 0.632, 0.815, 1), vec4(0.9, 1, 1, 1)},
 132      {vec4(0.999, 0.999, 0.519, 1), vec4(0.981, 0.667, 0.118, 1)}}};
 133 vec3 camPos = {0, 0, -3}, viewPos = {0, 0, 0};
 134 mat4 proj; // projection matrix
 135 mat4 view; // view matrix
 136 float znear = 0.001, zfar = 1000;
 137 float width = 512, height = 512;
 138 float aspect = width / height;
 139
 140 void setProjectionAndViewUniforms(GLuint progId) {
 141   GLuint projLoc = glGetUniformLocation(progId, "projection");
 142   glUniformMatrix4fv(projLoc, 1, GL_FALSE, glm::value_ptr(proj));
 143
 144   GLuint viewLoc = glGetUniformLocation(progId, "view");
 145   glUniformMatrix4fv(viewLoc, 1, GL_FALSE, glm::value_ptr(view));
 146 }
 147
 148 /** Orientates the transformation matrix to face the camera in the view matrix
 149  */
 150 mat4 faceView(mat4 m) {
 151   m[0][0] = view[0][0];
 152   m[0][1] = view[1][0];
 153   m[0][2] = view[2][0];
 154   m[1][0] = view[0][1];
 155   m[1][1] = view[1][1];
 156   m[1][2] = view[2][1];
 157   m[2][0] = view[0][2];
 158   m[2][1] = view[1][2];
 159   m[2][2] = view[2][2];
 160   return m;
 161 }
 162
 163 #ifdef PBO
 164 const int numPbos = 64;
 165 GLuint pboBufs[numPbos];
 166 GLbyte sink[1000 * 1000 * 4];
 167 #endif
 168
 169 void shadeClouds() {
 170   glDisable(GL_DEPTH_TEST);
 171   // shaderOutput * 0 + buffer * shader alpha
 172   glBlendFunc(GL_ZERO, GL_SRC_ALPHA);
 173   glEnable(GL_BLEND);
 174
 175   // sort by ascending distance from the sun
 176   sort(metaballs.begin(), metaballs.end(), [](Metaball &a, Metaball &b) {
 177     return distance(sunPos, a.pos) < distance(sunPos, b.pos);
 178   });
 179
 180   glActiveTexture(GL_TEXTURE0);
 181   glUniform1i(glGetUniformLocation(bbProg, "tex"), 0);
 182
 183   GLuint modelLoc = glGetUniformLocation(bbProg, "model");
 184   glUniform1i(glGetUniformLocation(bbProg, "debug"), 0);
 185
 186   /* glViewport(0, 0, shadeWidth, shadeHeight); */
 187
 188 #ifdef PBO
 189   int pboIdx = 0;
 190 #endif
 191
 192   auto begin_time = std::chrono::system_clock::now();
 193   size_t i = 0;
 194   for (auto &k : metaballs) {
 195     fprintf(stderr, "\rShading metaball %lu/%lu...", i++, metaballs.size());
 196     // place the billboard at the center of k
 197     mat4 model = translate(mat4(1), k.pos);
 198
 199     // rotate the billboard so that its normal is oriented to the sun
 200     model = faceView(model);
 201
 202     model = scale(model, vec3(metaballScale));
 203
 204     glUniformMatrix4fv(modelLoc, 1, GL_FALSE, glm::value_ptr(model));
 205
 206     // Set the billboard color as RGBA = (1.0, 1.0, 1.0, 1.0).
 207     vec4 color = {1, 1, 1, 1};
 208     glUniform4fv(glGetUniformLocation(bbProg, "color"), 1,
 209                  glm::value_ptr(color));
 210
 211     // Map the billboard texture with GL_MODULATE.
 212     // i.e. multiply rather than add
 213     // but glTexEnv is for the old fixed function pipeline --
 214     // need to just tell our fragment shader then to modulate
 215     int dIdx = k.d * NQ;
 216     glBindTexture(GL_TEXTURE_2D, bbTexIds[dIdx]);
 217     glUniform1i(glGetUniformLocation(bbProg, "modulate"), 1);
 218
 219     // Render the billboard.
 220     glDrawElements(GL_TRIANGLES, 6, GL_UNSIGNED_INT, 0);
 221
 222     // Read the pixel value corresponding to the center of metaball k.
 223     // 1. First get position in opengl screen space: from [-1,1]
 224     // 2. Normalize to [0,1]
 225     // 3. Multiply by (width * height)
 226     ivec2 screenPos =
 227         ((vec2(proj * view * model * vec4(0, 0, 0, 1)) + vec2(1)) / vec2(2)) *
 228         vec2(width, height);
 229     vec4 pixel;
 230
 231 #ifndef PBO
 232     // TODO: This is a huge bottleneck
 233     glReadPixels(screenPos.x, screenPos.y, 1, 1, GL_RGBA, GL_FLOAT, &pixel);
 234 #else
 235
 236     glBindBuffer(GL_PIXEL_PACK_BUFFER, pboBufs[pboIdx]);
 237     glReadPixels(0, 0, width, height, GL_RGBA, GL_UNSIGNED_BYTE,
 238                  NULL);
 239
 240     glBindBuffer(GL_PIXEL_PACK_BUFFER, pboBufs[(pboIdx + 1) % numPbos]);
 241     GLuint offset = screenPos.x * 4 + screenPos.y * (int)width * 4;
 242     ivec4 *src = (ivec4 *)glMapBufferRange(GL_PIXEL_PACK_BUFFER, offset, 4, GL_MAP_READ_BIT);
 243     checkError();
 244     if (src) {
 245       ivec4 t = src[(int)screenPos.x + (int)screenPos.y * (int)width];
 246       pixel = vec4(t) / vec4(255.f);
 247     }
 248     glUnmapBuffer(GL_PIXEL_PACK_BUFFER);
 249
 250     pboIdx = (pboIdx + 1) % numPbos;
 251
 252     glBindBuffer(GL_PIXEL_PACK_BUFFER, 0);
 253 #endif
 254
 255     // Multiply the pixel value by the sunlight color.
 256     pixel *= envColors[envColorIdx][0];
 257
 258     // Store the color into an array C[k] as the color of the billboard.
 259     k.col = pixel;
 260   }
 261   fprintf(stderr, "\n");
 262
 263   auto elapsed = std::chrono::system_clock::now() - begin_time;
 264   double elapsed_seconds = std::chrono::duration_cast<std::chrono::duration<double> >(elapsed).count();
 265   fprintf(stderr, "time with pbo: %f\n", elapsed_seconds);
 266
 267   saveFBO();
 268   checkError();
 269   /* glViewport(0, 0, width, height); */
 270 }
 271
 272 void renderObject() {
 273   glDisable(GL_BLEND);
 274   // render the sun
 275   glUseProgram(sunProg);
 276   mat4 model = translate(mat4(1), sunPos);
 277   /* model = lookAt(sunPos, sunPos + sunDir, {0, 1, 0}) * model; */
 278   model = translate(scale(translate(model, -sunPos), vec3(0.3)), sunPos);
 279   glUniformMatrix4fv(glGetUniformLocation(sunProg, "model"), 1, GL_FALSE,
 280                      glm::value_ptr(model));
 281   glDrawElements(GL_TRIANGLES, 6, GL_UNSIGNED_INT, 0);
 282 }
 283
 284 void renderClouds() {
 285   glUseProgram(bbProg);
 286
 287   // Sort metaballs in descending order from the viewpoint
 288   sort(metaballs.begin(), metaballs.end(), [](Metaball &a, Metaball &b) {
 289     return distance(camPos, a.pos) > distance(camPos, b.pos);
 290   });
 291
 292   glUniform1i(glGetUniformLocation(bbProg, "debug"), curMode != render);
 293
 294   glDisable(GL_DEPTH_TEST);
 295   glEnable(GL_BLEND);
 296   // shaderOutput * 1 + buffer * shader alpha
 297   glBlendFunc(GL_ONE, GL_SRC_ALPHA);
 298
 299   /* glBlendColor(1.f,1.f,1.f,1.f); */
 300   /* glBlendFuncSeparate(GL_ONE, GL_SRC_ALPHA, GL_CONSTANT_ALPHA, GL_SRC_ALPHA);
 301    */
 302
 303   glActiveTexture(GL_TEXTURE0);
 304   glUniform1i(glGetUniformLocation(bbProg, "tex"), 0);
 305
 306   for (int i = 0; i < metaballs.size(); i++) {
 307     Metaball k = metaballs[i];
 308
 309     GLuint modelLoc = glGetUniformLocation(bbProg, "model");
 310
 311     // Place the billboard at the center of the corresponding metaball n.
 312     mat4 model = translate(mat4(1), k.pos);
 313     // Rotate the billboard so that its normal is oriented to the viewpoint.
 314     model = faceView(model);
 315
 316     model = scale(model, vec3(metaballScale));
 317
 318     glUniformMatrix4fv(modelLoc, 1, GL_FALSE, glm::value_ptr(model));
 319
 320     // Set the billboard color as C[n].
 321     k.col.w = 1;
 322     glUniform4fv(glGetUniformLocation(bbProg, "color"), 1,
 323                  glm::value_ptr(k.col));
 324
 325     // Map the billboard texture.
 326     int dIdx = k.d * (NQ - 1);
 327     glBindTexture(GL_TEXTURE_2D, bbTexIds[dIdx]);
 328
 329     // Don't modulate it -- blend it
 330     glUniform1i(glGetUniformLocation(bbProg, "modulate"), 0);
 331
 332     glUniform1f(glGetUniformLocation(bbProg, "debugColor"),
 333                 curMode == debugColor);
 334     if (curMode != render) {
 335       float debugVal = 0;
 336       if (curMode == debugContDist)
 337         debugVal = k.d;
 338       else if (curMode == debugProbAct)
 339         debugVal = cs.p_act[k.coords.x][k.coords.y][k.coords.z] / P_ACT;
 340       else if (curMode == debugProbExt)
 341         debugVal = cs.p_ext[k.coords.x][k.coords.y][k.coords.z] / P_EXT;
 342       glUniform1f(glGetUniformLocation(bbProg, "debugVal"), debugVal);
 343       glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
 344       model = scale(model, vec3(0.2));
 345       glUniformMatrix4fv(modelLoc, 1, GL_FALSE, glm::value_ptr(model));
 346     }
 347
 348     // Render the billboard with the blending function.
 349     glDrawElements(GL_TRIANGLES, 6, GL_UNSIGNED_INT, 0);
 350   }
 351 }
 352
 353 bool needsReshading = true;
 354 void display() {
 355   if (needsReshading) {
 356     // TODO: find a way to make sure there's no clipping
 357     view = glm::lookAt(sunPos + sunDir * vec3(100.f), sunPos, {0, 0, 1});
 358     // TODO: calculate bounds so everything is covered
 359     proj = glm::ortho(2.5f, -2.5f, -2.5f, 2.5f, znear, 10000.f);
 360     glUseProgram(bbProg);
 361     setProjectionAndViewUniforms(bbProg);
 362
 363     glClearColor(1, 1, 1, 1);
 364     glClear(GL_COLOR_BUFFER_BIT);
 365     shadeClouds();
 366     needsReshading = false;
 367   }
 368
 369   view = glm::lookAt(camPos, viewPos, {0, 1, 0});
 370   proj = glm::perspective(45.f, aspect, znear, zfar);
 371   glUseProgram(sunProg);
 372   setProjectionAndViewUniforms(sunProg);
 373   glUseProgram(bbProg);
 374   setProjectionAndViewUniforms(bbProg);
 375
 376   vec4 skyColor = envColors[envColorIdx][1];
 377   glClearColor(skyColor.r, skyColor.g, skyColor.b,
 378                skyColor.a); // background color
 379   glClear(GL_COLOR_BUFFER_BIT);
 380   renderObject(); // render things that aren't clouds
 381   renderClouds();
 382
 383   glutSwapBuffers();
 384 }
 385
 386 bool needsRedisplay = false;
 387 void timer(int _) {
 388   /* calculateMetaballs(); */
 389   if (needsRedisplay) {
 390     glutPostRedisplay();
 391   }
 392   needsRedisplay = false;
 393   glutTimerFunc(16, timer, 0);
 394 }
 395
 396 void keyboard(unsigned char key, int x, int y) {
 397   if (key == ' ') {
 398     calculateMetaballs();
 399     needsRedisplay = true;
 400     needsReshading = curMode == render;
 401   }
 402   if (key == '0') {
 403     needsReshading = curMode != render;
 404     curMode = render;
 405     needsRedisplay = true;
 406   }
 407   if (key == '1') {
 408     curMode = debugContDist;
 409     needsRedisplay = true;
 410   }
 411   if (key == '2') {
 412     curMode = debugColor;
 413     needsRedisplay = true;
 414   }
 415   if (key == '3') {
 416     curMode = debugProbAct;
 417     needsRedisplay = true;
 418   }
 419   if (key == '4') {
 420     curMode = debugProbExt;
 421     needsRedisplay = true;
 422   }
 423   if (key == 's') {
 424     envColorIdx = (envColorIdx + 1) % envColors.size();
 425     needsRedisplay = true;
 426     needsReshading = true;
 427   }
 428 }
 429
 430 int prevMouseX, prevMouseY;
 431 bool firstMouse = true;
 432 void motion(int x, int y) {
 433   if (firstMouse) {
 434     prevMouseX = x;
 435     prevMouseY = y;
 436     firstMouse = false;
 437   }
 438   float dx = x - prevMouseX, dy = y - prevMouseY;
 439   prevMouseX = x;
 440   prevMouseY = y;
 441   const vec3 origin(0, 0, 0);
 442   const float sensitivity = 0.003f;
 443   auto camMat = translate(mat4(1), origin + camPos);
 444   auto rotation = rotate(rotate(mat4(1), -dx * sensitivity, {0, 1, 0}),
 445                          -dy * sensitivity, {1, 0, 0});
 446   auto rotAroundOrig = camMat * rotation * translate(mat4(1), origin - camPos);
 447   camPos = rotAroundOrig * glm::vec4(camPos, 0);
 448   needsRedisplay = true;
 449 }
 450
 451 void passiveMotion(int x, int y) {
 452   prevMouseX = x;
 453   prevMouseY = y;
 454 }
 455
 456 int main(int argc, char **argv) {
 457   glutInit(&argc, argv);
 458   glutInitDisplayMode(GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGBA | GLUT_ALPHA |
 459                       GLUT_3_2_CORE_PROFILE);
 460   glutInitWindowSize(width, height);
 461   glutCreateWindow("Clouds");
 462   glutDisplayFunc(display);
 463
 464   glewInit();
 465
 466   Program prog("billboardvert.glsl", "billboardfrag.glsl");
 467   bbProg = prog.progId;
 468   Program sProg("sunvert.glsl", "sunfrag.glsl");
 469   sunProg = sProg.progId;
 470
 471   glGenVertexArrays(1, &bbVao);
 472   glUseProgram(sunProg);
 473   glBindVertexArray(bbVao);
 474   glUseProgram(bbProg);
 475   glBindVertexArray(bbVao);
 476   GLuint vbos[2];
 477   glGenBuffers(2, vbos);
 478
 479   vector<vec3> poss = {{-1, -1, 0}, {-1, 1, 0}, {1, 1, 0}, {1, -1, 0}};
 480   vector<GLuint> indices = {2, 1, 0, 3, 2, 0};
 481
 482   GLuint posLoc = glGetAttribLocation(bbProg, "vPosition");
 483   glBindBuffer(GL_ARRAY_BUFFER, vbos[0]);
 484   glBufferData(GL_ARRAY_BUFFER, poss.size() * sizeof(glm::vec3), &poss[0],
 485                GL_STATIC_DRAW);
 486   glEnableVertexAttribArray(posLoc);
 487   glVertexAttribPointer(posLoc, 3, GL_FLOAT, GL_FALSE, 0, 0);
 488
 489   glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, vbos[1]);
 490   glBufferData(GL_ELEMENT_ARRAY_BUFFER, indices.size() * sizeof(GLuint),
 491                &indices[0], GL_STATIC_DRAW);
 492
 493   prog.validate();
 494
 495   precalculateBillboardTextures();
 496
 497   initClouds(&cs);
 498   calculateMetaballs();
 499
 500   glutKeyboardFunc(keyboard);
 501   glutMotionFunc(motion);
 502   glutPassiveMotionFunc(passiveMotion);
 503   glutTimerFunc(16, timer, 0);
 504
 505 #ifdef PBO
 506   // setup PBOs for buffering readPixels
 507   glGenBuffers(numPbos, pboBufs);
 508   for (int i = 0; i < numPbos; i++) {
 509     glBindBuffer(GL_PIXEL_PACK_BUFFER, pboBufs[i]);
 510     glBufferData(GL_PIXEL_PACK_BUFFER, width * height * 4, NULL,
 511                  GL_DYNAMIC_READ);
 512   }
 513   glBindBuffer(GL_PIXEL_PACK_BUFFER, 0);
 514 #endif
 515
 516   glutMainLoop();
 517
 518   return 0;
 519 }