Continuous distribution
authorLuke Lau <luke_lau@icloud.com>
Fri, 10 Apr 2020 14:02:27 +0000 (15:02 +0100)
committerLuke Lau <luke_lau@icloud.com>
Fri, 10 Apr 2020 14:02:27 +0000 (15:02 +0100)
.gitignore
Makefile
billboardfrag.glsl
clouds.cpp
debug.hpp
simulation.cpp
simulation.h

index dfe6a70738615a3dd09c385f29e3639b1e2203b1..309bd5cda988895176ce6a33aff1bb3284d9ea3b 100644 (file)
@@ -3,3 +3,4 @@
 clouds
 fbo.tga
 *.dSYM
+bbtex
index 67b688d8f657a20bf428114d10798df30a3f8cb8..64658142a5a8bae7bf4c647ccf3c5b57a9f59cc3 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,5 +1,5 @@
 clouds: clouds.cpp program.cpp simulation.cpp
-       clang++ -std=c++17 $^ -o $@ \
+       clang++ -std=c++17 $^ -march=native -o $@ \
                -Wall -g \
                -I/usr/local/include -L/usr/local/lib \
                -framework OpenGL -framework glut -lglew
index 94d15dcb2f1e477e90ed7cba724dedfcc4f49308..046b7bd3bca6e335cd847e6e6dc4151f4edf03fa 100644 (file)
@@ -4,7 +4,13 @@ uniform sampler2D tex;
 in vec2 texCoord;
 out vec4 FragColor;
 uniform bool modulate;
+uniform bool debugContDist;
+uniform float density;
 void main() {
+  if (debugContDist) {
+    FragColor = mix(vec4(1, 1, 1, 1), vec4(1, 0, 0, 1), density);
+    return;
+  }
        // Cf = color from fragment, Ct = color from texture
        // Cc = color from texture environment -- not set, defaults to (0,0,0,0)
        // Af = alpha from fragment, At = alpha from texture
index d5c8f6ad628ff1e1c23d22a09abfdd5e7a180bcd..52b3053812d8f16416361e519270e8df1bd23dfa 100644 (file)
@@ -1,27 +1,26 @@
 #include "debug.hpp"
-#include "simulation.h"
 #include "program.hpp"
+#include "simulation.h"
 #include <GL/glew.h>
 #include <GLUT/glut.h>
+#include <array>
 #include <cmath>
 #include <cstdio>
 #include <cstdlib>
 #include <glm/ext.hpp>
 #include <glm/glm.hpp>
+#include <sys/stat.h>
 #include <vector>
 
 #pragma clang diagnostic ignored "-Wdeprecated-declarations"
 
+bool debugContDist = false;
+
 using namespace std;
 using namespace glm;
 
-// calculate continuous distribution
-void calcContDist(Clouds *clds);
-
-float w(int i, int j, int k) { return 1; }
-
-const float metaballR = 1;
-float metaballField(float r) {
+const float metaballR = 1.5f;
+inline float metaballField(float r) {
   if (r > metaballR)
     return 0;
   const float a = r / metaballR;
@@ -29,34 +28,7 @@ float metaballField(float r) {
          (22.f / 9.f * powf(a, 2)) + 1;
 }
 
-/* const float normalizationFactor = 748.f / 405.f * M_PI * metaballR; */
-
-void calcContDist(Clouds *clds, float t) {
-  const int i0 = 2, j0 = 2, k0 = 2, t0 = 2;
-  const float divisor =
-      1.f / ((2 * t0 + 1) * (2 * k0 + 1) * (2 * j0 + 1) * (2 * i0 + 1));
-  float sum = 0;
-  for (int i = 0; i < CLOUD_DIM; i++) {
-    for (int j = 0; j < CLOUD_DIM; j++) {
-      for (int k = 0; k < CLOUD_DIM; k++) {
-
-        // inner sums
-        /* for (int tp = -t0, tp < t0; tp++) { */
-        for (int ip = -i0; ip < i0; ip++) {
-          for (int jp = -j0; jp < j0; jp++) {
-            for (int kp = -k0; kp < k0; kp++) {
-
-              sum += w(ip, jp, kp) * (float)clds->cld[i + ip][j + jp][k + kp];
-            }
-          }
-        }
-        /* } */
-
-        clds->contDist[i][j][k] = sum / divisor;
-      }
-    }
-  }
-}
+const float normalizationFactor = 748.f / 405.f * M_PI * metaballR;
 
 void checkError() {
   if (GLenum e = glGetError()) {
@@ -65,8 +37,6 @@ void checkError() {
   }
 }
 
-vector<vec4> bbColors;
-
 GLuint bbProg;
 GLuint bbVao;
 
@@ -74,67 +44,78 @@ GLuint bbVao;
 // These textures then go on the billboards
 // The texture stores attenuation ratio?
 
-#define NQ 1
+#define NQ 64
 GLuint bbTexIds[NQ];
 
 // Stores attenuation ratio inside r channel
 // Should be highest value at center
 void precalculateBillboardTextures() {
+  fprintf(stderr, "Calculating billboard textures...\n");
+  glGenTextures(NQ, bbTexIds);
+
+  for (int d = 0; d < NQ; d++) {
     float data[32 * 32];
+    for (int j = 0; j < 32; j++) {
+      for (int i = 0; i < 32; i++) {
         // TODO: properly calculate this instead of whatever this is
-  for (int j = 0; j < 32; j++)
-    for (int i = 0; i < 32; i++)
-      data[i + j * 32] = fmin(1.f, 0.5f + 2.f * (distance(vec2(i, j), vec2(16, 16)) / 16));
+        float r = distance(vec2(i, j), vec2(16, 16)) / 16;
+        float density = (float)d / NQ;
+        data[i + j * 32] =
+            1 - (density * metaballField(r) / normalizationFactor);
+      }
+    }
 
-  glGenTextures(NQ, bbTexIds);
+    mkdir("bbtex", 0777);
+    char path[32];
+    snprintf(path, 32, "bbtex/%i.tga", d);
+    saveGrayscale(data, 32, 32, path);
 
-  for (int i = 0; i < NQ; i++) {
-    glBindTexture(GL_TEXTURE_2D, bbTexIds[i]);
+    glBindTexture(GL_TEXTURE_2D, bbTexIds[d]);
     checkError();
 
     glTexImage2D(GL_TEXTURE_2D, 0, GL_RED, 32, 32, 0, GL_RED, GL_FLOAT, data);
     glGenerateMipmap(GL_TEXTURE_2D); // required, otherwise texture is blank
 
     checkError();
+
+    fprintf(stderr, "\r%i out of %i densities calculated%s", d + 1, NQ,
+            d == NQ - 1 ? "\n" : "");
   }
 }
 
 struct Metaball {
   vec3 pos;
-  float r;
+  /** Radius, density */
+  float r, d;
+  vec4 col;
 };
-// TODO: why is the x axis flipped??
-/* vector<Metaball> metaballs = {{{-0.5, 0.5, 0.5}, 0.25}, */
-/*                               {{-0.3, 0.5, 0.3}, 0.25}}; */
-vector<Metaball> metaballs = {{{0, 0, 0.5}, 1.f},
-                              {{0, 0.3, 0.3}, 0.7f}};
+
+array<Metaball, CLOUD_DIM * CLOUD_DIM * CLOUD_DIM> metaballs;
 
 Clouds cs;
 
 void calculateMetaballs() {
-  /* stepClouds(&cs); */
-  metaballs.clear();
-  for (int i = 0; i < 256; i++) {
-    float x = ((float)rand()/(float)(RAND_MAX) - 0.5) * 2;
-    float y = ((float)rand()/(float)(RAND_MAX) - 0.5) * 2;
-    float z = ((float)rand()/(float)(RAND_MAX) - 0.5) * 2;
-    float r = (float)rand()/(float)(RAND_MAX) * 1;
-    Metaball m = {{x,y,z}, r};
-    metaballs.push_back(m);
-  }
-  /* for (int i = 0; i < CLOUD_DIM; i++) { */
-  /*   for (int j = 0; j < CLOUD_DIM; j++) { */
-  /*     for (int k = 0; k < CLOUD_DIM; k++) { */
-  /*       if (cs.cld[i][j][k]) { */
-  /*         Metaball m = {{i / (float)CLOUD_DIM, j / (float)CLOUD_DIM, k / (float)CLOUD_DIM}, */
-  /*                       1.f / (float)CLOUD_DIM }; */
-  /*         m.pos = (m.pos * vec3(2)) - vec3(1); */
+  stepClouds(&cs);
+  /* for (int i = 0; i < 256; i++) { */
+  /*   float x = ((float)rand()/(float)(RAND_MAX) - 0.5) * 2; */
+  /*   float y = ((float)rand()/(float)(RAND_MAX) - 0.5) * 2; */
+  /*   float z = ((float)rand()/(float)(RAND_MAX) - 0.5) * 2; */
+  /*   float r = (float)rand()/(float)(RAND_MAX) * 1; */
+  /*   Metaball m = {{x,y,z}, r}; */
   /*   metaballs.push_back(m); */
   /* } */
-  /*     } */
-  /*   } */
-  /* } */
-  fprintf(stderr, "num metaballs: %lu\n", metaballs.size());
+  for (int i = 0; i < CLOUD_DIM; i++) {
+    for (int j = 0; j < CLOUD_DIM; j++) {
+      for (int k = 0; k < CLOUD_DIM; k++) {
+        Metaball m = {
+            {i / (float)CLOUD_DIM, j / (float)CLOUD_DIM, k / (float)CLOUD_DIM},
+            1.f / (float)CLOUD_DIM};
+        m.pos = (m.pos * vec3(2)) - vec3(1);
+        m.d = cs.q[i][j][k];
+        metaballs[i * CLOUD_DIM * CLOUD_DIM + j * CLOUD_DIM + k] = m;
+      }
+    }
+  }
 }
 
 vec3 sunPos = {0, 2, 2}, sunDir = {0, -1, -1};
@@ -171,9 +152,6 @@ mat4 faceView(mat4 m) {
 GLuint attenuationTex;
 
 void shadeClouds() {
-
-  bbColors.clear();
-
   glDisable(GL_DEPTH_TEST);
   // shaderOutput * 0 + buffer * shader alpha
   glBlendFunc(GL_ZERO, GL_SRC_ALPHA);
@@ -185,12 +163,11 @@ void shadeClouds() {
   });
 
   glActiveTexture(GL_TEXTURE0);
-  glBindTexture(GL_TEXTURE_2D, bbTexIds[0]);
   glUniform1i(glGetUniformLocation(bbProg, "tex"), 0);
 
   GLuint modelLoc = glGetUniformLocation(bbProg, "model");
 
-  for (auto k : metaballs) {
+  for (auto &k : metaballs) {
     // place the billboard at the center of k
     mat4 model = scale(translate(mat4(1), k.pos), vec3(k.r) * 2.f);
 
@@ -208,6 +185,8 @@ void shadeClouds() {
     // i.e. multiply rather than add
     // but glTexEnv is for the old fixed function pipeline --
     // need to just tell our fragment shader then to modulate
+    int dIdx = k.d * NQ;
+    glBindTexture(GL_TEXTURE_2D, bbTexIds[dIdx]);
     glUniform1i(glGetUniformLocation(bbProg, "modulate"), 1);
 
     // Render the billboard.
@@ -217,20 +196,20 @@ void shadeClouds() {
     // 1. First get position in opengl screen space: from [-1,1]
     // 2. Normalize to [0,1]
     // 3. Multiply by (width * height)
-    vec2 screenPos = ((vec2(proj * view * model * vec4(0,0,0,1)) + vec2(1)) / vec2(2))
-                      * vec2(width, height);
+    vec2 screenPos =
+        ((vec2(proj * view * model * vec4(0, 0, 0, 1)) + vec2(1)) / vec2(2)) *
+        vec2(width, height);
     vec4 pixel;
-    glReadPixels(screenPos.x, screenPos.y, 1, 1, GL_RGBA, GL_FLOAT, value_ptr(pixel));
-    /* if (pixel.g == 0 && pixel.b == 0) abort(); */
-    /* fprintf(stderr, "pixel:"); */
-    /* dump(pixel); */
+    // TODO: This is a huge bottleneck
+    glReadPixels(screenPos.x, screenPos.y, 1, 1, GL_RGBA, GL_FLOAT,
+                 value_ptr(pixel));
 
     // Multiply the pixel value by the sunlight color.
     vec4 sunColor = {1, 1, 0.9, 1};
     pixel *= sunColor;
 
     // Store the color into an array C[k] as the color of the billboard.
-    bbColors.push_back(pixel);
+    k.col = pixel;
   }
 
   saveFBO();
@@ -245,10 +224,16 @@ void renderClouds() {
     return distance(camPos, a.pos) > distance(camPos, b.pos);
   });
 
+  glUniform1i(glGetUniformLocation(bbProg, "debugContDist"), debugContDist);
+
   glDisable(GL_DEPTH_TEST);
   glEnable(GL_BLEND);
   // shaderOutput * 1 + buffer * shader alpha
   glBlendFunc(GL_ONE, GL_SRC_ALPHA);
+
+  glActiveTexture(GL_TEXTURE0);
+  glUniform1i(glGetUniformLocation(bbProg, "tex"), 0);
+
   for (int i = 0; i < metaballs.size(); i++) {
     Metaball k = metaballs[i];
 
@@ -262,29 +247,32 @@ void renderClouds() {
     glUniformMatrix4fv(modelLoc, 1, GL_FALSE, glm::value_ptr(model));
 
     // Set the billboard color as C[n].
-    fprintf(stderr, "bbColors[i]: ");
-    dump(bbColors[i]);
-    /* bbColors[i].x = 1 - bbColors[i].x; */
-    /* bbColors[i].y = 1 - bbColors[i].y; */
-    /* bbColors[i].z = 1 - bbColors[i].z; */
-    bbColors[i].w = 1;
+    k.col.w = 1;
     glUniform4fv(glGetUniformLocation(bbProg, "color"), 1,
-                 glm::value_ptr(bbColors[i]));
+                 glm::value_ptr(k.col));
 
     // Map the billboard texture.
-    glActiveTexture(GL_TEXTURE0);
-    glBindTexture(GL_TEXTURE_2D, bbTexIds[0]);
-    glUniform1i(glGetUniformLocation(bbProg, "tex"), 0);
+    int dIdx = k.d * NQ;
+    glBindTexture(GL_TEXTURE_2D, bbTexIds[dIdx]);
 
     // Don't modulate it -- blend it
     glUniform1i(glGetUniformLocation(bbProg, "modulate"), 0);
 
+    if (debugContDist) {
+      glUniform1f(glGetUniformLocation(bbProg, "density"), k.d);
+      glDisable(GL_BLEND);
+      model = scale(model, vec3(0.02));
+      glUniformMatrix4fv(modelLoc, 1, GL_FALSE, glm::value_ptr(model));
+    }
+
     // Render the billboard with the blending function.
     glDrawElements(GL_TRIANGLES, 6, GL_UNSIGNED_INT, 0);
   }
 }
 
+bool needsReshading = true;
 void display() {
+  if (needsReshading) {
     // TODO: find a way to make sure there's no clipping
     view = glm::lookAt(sunPos + sunDir * vec3(20), sunPos, {0, 1, 0});
     proj = glm::ortho(1.f * aspect, -1.f * aspect, -1.f, 1.f, znear, zfar);
@@ -293,6 +281,8 @@ void display() {
     glClearColor(1, 1, 1, 1);
     glClear(GL_COLOR_BUFFER_BIT);
     shadeClouds();
+    needsReshading = false;
+  }
 
   view = glm::lookAt(camPos, viewPos, {0, 1, 0});
   proj = glm::perspective(45.f, aspect, znear, zfar);
@@ -319,7 +309,12 @@ void timer(int _) {
 void keyboard(unsigned char key, int x, int y) {
   if (key == ' ') {
     calculateMetaballs();
-               glutPostRedisplay();
+    needsRedisplay = true;
+    needsReshading = true;
+  }
+  if (key == 'd') {
+    debugContDist = !debugContDist;
+    needsRedisplay = true;
   }
 }
 
@@ -332,7 +327,8 @@ void motion(int x, int y) {
     firstMouse = false;
   }
   float dx = x - prevMouseX, dy = y - prevMouseY;
-       prevMouseX = x; prevMouseY = y;
+  prevMouseX = x;
+  prevMouseY = y;
   const vec3 origin(0, 18, 0);
   const float sensitivity = 0.003f;
   auto camMat = translate(mat4(1), origin + camPos);
@@ -343,6 +339,11 @@ void motion(int x, int y) {
   needsRedisplay = true;
 }
 
+void passiveMotion(int x, int y) {
+  prevMouseX = x;
+  prevMouseY = y;
+}
+
 int main(int argc, char **argv) {
   glutInit(&argc, argv);
   glutInitDisplayMode(GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGB |
@@ -388,6 +389,7 @@ int main(int argc, char **argv) {
 
   glutKeyboardFunc(keyboard);
   glutMotionFunc(motion);
+  glutPassiveMotionFunc(passiveMotion);
   glutTimerFunc(16, timer, 0);
 
   // set up billboard prog
index debe65e6eab44875f8ed97e02b8f837ae6a8aeaf..6774a1a66c6724fc0e9d35c0516e23dd28e2c5e5 100644 (file)
--- a/debug.hpp
+++ b/debug.hpp
@@ -20,7 +20,7 @@ inline void dump(glm::vec2 v) { fprintf(stderr, "{%f,%f}\n", v.x, v.y); }
 typedef struct {
   char head[12];
   short dx /* Width */, dy /* Height */, head2;
-  char pic[600 * 400][3];
+  char pic[600 * 400 * 3];
 } TGA;
 char tgahead[12] = {0x00, 0x00, 0x02, 0x00, 0x00, 0x00,
                     0x00, 0x00, 0x00, 0x00, 0x00, 0x00};
@@ -33,9 +33,36 @@ void saveFBO() {
   tga->dy = height;
   tga->head2 = 0x2018;
   // TODO: flip -- the image is upside down so far
-  glReadPixels(0, 0, width, height, GL_BGR, GL_UNSIGNED_BYTE, tga->pic[0]);
+  glReadPixels(0, 0, width, height, GL_BGR, GL_UNSIGNED_BYTE, tga->pic);
   FILE *cc = fopen("fbo.tga", "wb");
   fwrite(tga, 1, (18 + 3 * width * height), cc);
   fclose(cc);
   free(tga);
 }
+
+void saveImage(char *img, short width, short height, const char *file) {
+  TGA *tga = (TGA *)malloc(sizeof(TGA));
+  memcpy(tga->head, tgahead, 12);
+  tga->dx = width;
+  tga->dy = height;
+  tga->head2 = 0x2018;
+  memcpy(tga->pic, img, width * height * 3);
+  // TODO: flip -- the image is upside down so far
+  FILE *cc = fopen(file, "wb");
+  fwrite(tga, 1, (18 + 3 * width * height), cc);
+  fclose(cc);
+  free(tga);
+}
+
+void saveGrayscale(float *img, short width, short height, const char *file) {
+  char *bytes = (char*)malloc(sizeof(char) * width * height * 3);
+  for (int i = 0; i < width; i++) {
+    for (int j = 0; j < height; j++) {
+      for (int z = 0; z < 3; z++) {
+        bytes[i * 3 + j * width * 3 + z] = (int)(img[i + j * width] * 255.f * 2);
+      }
+    }
+  }
+  saveImage(bytes, width, height, file);
+  free(bytes);
+}
index 0471f4a20b4cba1dba112237513f940df07051fe..e022d12b1b921901b42dcceaa50466a0e96b439f 100644 (file)
@@ -98,7 +98,43 @@ void extinction(Clouds *cs) {
   *cs = ncs;
 }
 
+/** Weighting function */
+float w(int ip, int jp, int kp) { return 1; }
+
+void calcContDist(Clouds *cls) {
+  const int i0 = 2, j0 = 2, k0 = 2, t0 = 0;
+  const float divisor =
+      1.f / ((2 * t0 + 1) * (2 * k0 + 1) * (2 * j0 + 1) * (2 * i0 + 1));
+  for (int i = 0; i < CLOUD_DIM; i++) {
+    for (int j = 0; j < CLOUD_DIM; j++) {
+      for (int k = 0; k < CLOUD_DIM; k++) {
+        float sum = 0;
+
+        // sum
+        for (int tp = -t0; tp <= t0; tp++) {
+          for (int ip = -i0; ip <= i0; ip++) {
+            for (int jp = -j0; jp <= j0; jp++) {
+              for (int kp = -k0; kp <= k0; kp++) {
+                if (i + ip < 0 || i + ip >= CLOUD_DIM ||
+                    j + jp < 0 || j + jp >= CLOUD_DIM ||
+                    k + kp < 0 || k + kp >= CLOUD_DIM)
+                  continue;
+
+                sum += w(ip, jp, kp) * (float)cls->cld[i + ip][j + jp][k + kp];
+              }
+            }
+          }
+        }
+
+        cls->q[i][j][k] = sum * divisor;
+        if (cls->q[i][j][k] > 1) abort();
+      }
+    }
+  }
+}
+
 void stepClouds(Clouds *cs) {
   growth(cs);
   extinction(cs);
+  calcContDist(cs);
 }
index f021ce4b6c31f58f389af55381d08efb7c915be7..c6c9e03b02cae0ae2938124685750bf7e8a793da 100644 (file)
@@ -7,7 +7,8 @@ struct Clouds {
   float p_ext[CLOUD_DIM][CLOUD_DIM][CLOUD_DIM];
   float p_hum[CLOUD_DIM][CLOUD_DIM][CLOUD_DIM];
   float p_act[CLOUD_DIM][CLOUD_DIM][CLOUD_DIM];
-  float contDist[CLOUD_DIM][CLOUD_DIM][CLOUD_DIM];
+  /** continuous distribution */
+  float q[CLOUD_DIM][CLOUD_DIM][CLOUD_DIM];
 };
 
 void initClouds(Clouds *cs);