clouds.cpp

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