Draw the sun, debug the colour and variable size
authorLuke Lau <luke_lau@icloud.com>
Fri, 10 Apr 2020 23:19:12 +0000 (00:19 +0100)
committerLuke Lau <luke_lau@icloud.com>
Fri, 10 Apr 2020 23:19:12 +0000 (00:19 +0100)
billboardfrag.glsl
clouds.cpp
simulation.cpp
simulation.h
sunfrag.glsl [new file with mode: 0644]
sunvert.glsl [new file with mode: 0644]

index 16e074e616aae9ea8101b38a11564982e249c894..c2e480bf88a930278ffe6f472f8be86c02b1bddd 100644 (file)
@@ -5,8 +5,13 @@ in vec2 texCoord;
 out vec4 FragColor;
 uniform bool modulate;
 uniform bool debug;
+uniform bool debugColor;
 uniform float debugVal;
 void main() {
+  if (debugColor) {
+    FragColor = color;
+    return;
+  }
   if (debug) {
     FragColor = mix(vec4(1, 1, 1, 1), vec4(1, 0, 0, 1), debugVal);
     return;
@@ -22,6 +27,9 @@ void main() {
                // A = Af * At                          
                FragColor = color * f;
        } else {
+        // "That is, the colors in the frame buffer are multiplied by the
+        // attenuation ratio of the billboard texture and then the colors in
+        // the texture are added"
                // GL_BLEND:
                // C = Cf * (1-Ct) + Cc * Ct
                // A = Af * At
index baabcba703d4aca0e44e80f142e7eb97a82ed1c7..5d188233fe66c9eb16bdbbfcd5a08fffd7a377cb 100644 (file)
@@ -17,6 +17,7 @@
 enum Mode {
   render,
   debugContDist,
+  debugColor,
   debugProbExt,
   debugProbAct
 };
@@ -25,7 +26,7 @@ Mode curMode = render;
 using namespace std;
 using namespace glm;
 
-const float metaballR = 1.5f;
+const float metaballR = 2.f * 1.f / 16.f;
 inline float metaballField(float r) {
   if (r > metaballR)
     return 0;
@@ -43,7 +44,7 @@ void checkError() {
   }
 }
 
-GLuint bbProg;
+GLuint bbProg, sunProg;
 GLuint bbVao;
 
 // Here we need to generate n_q textures for different densities of metaballs
@@ -64,10 +65,10 @@ void precalculateBillboardTextures() {
     for (int j = 0; j < 32; j++) {
       for (int i = 0; i < 32; i++) {
         // TODO: properly calculate this instead of whatever this is
-        float r = distance(vec2(i, j), vec2(16, 16)) / 16;
+        float r = distance(vec2(i, j), vec2(16, 16)) / 32;
         float density = (float)d / NQ;
         data[i + j * 32] =
-            1 - (density * metaballField(r) / normalizationFactor);
+            1 - (density * 0.01 * metaballField(r * metaballR) / normalizationFactor);
       }
     }
 
@@ -92,12 +93,12 @@ void precalculateBillboardTextures() {
 struct Metaball {
   vec3 pos;
   ivec3 coords;
-  /** Radius, density */
-  float r, d;
+  /** Density */
+  float d;
   vec4 col;
 };
 
-array<Metaball, CLOUD_DIM * CLOUD_DIM * CLOUD_DIM> metaballs;
+array<Metaball, CLOUD_DIM_X * CLOUD_DIM_Y * CLOUD_DIM_Z> metaballs;
 
 Clouds cs;
 
@@ -111,22 +112,25 @@ void calculateMetaballs() {
   /*   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++) {
+  for (int i = 0; i < CLOUD_DIM_X; i++) {
+    for (int j = 0; j < CLOUD_DIM_Y; j++) {
+      for (int k = 0; k < CLOUD_DIM_Z; k++) {
+        const float cloudScale = 1.f / 16;
         Metaball m = {
-            {i / (float)CLOUD_DIM, j / (float)CLOUD_DIM, k / (float)CLOUD_DIM},
-            {i, j, k},
-            1.f / (float)CLOUD_DIM};
-        m.pos = (m.pos * vec3(2)) - vec3(1);
+            vec3(i, j, k) * vec3(cloudScale),
+            {i, j, k} };
+        /* m.pos = (m.pos * vec3(2)) - (cloudScale / 2); */
+        m.pos -= vec3(CLOUD_DIM_X, CLOUD_DIM_Y, CLOUD_DIM_Z) * cloudScale / 2.f;
         m.d = cs.q[i][j][k];
-        metaballs[i * CLOUD_DIM * CLOUD_DIM + j * CLOUD_DIM + k] = m;
+        metaballs[i * CLOUD_DIM_Y * CLOUD_DIM_Z + j * CLOUD_DIM_Z + k] = m;
       }
     }
   }
 }
 
-vec3 sunPos = {0, 2, 2}, sunDir = {0, -1, -1};
+vec3 sunPos = {0, 5, 5}, sunDir = {0, -1, -1};
+/* vec4 sunColor = {1,0,0.429,1}; */
+vec4 sunColor = {1,1,1,1};
 vec3 camPos = {0, 0, -5}, viewPos = {0, 0, 0};
 mat4 proj; // projection matrix
 mat4 view; // view matrix
@@ -159,6 +163,8 @@ mat4 faceView(mat4 m) {
 
 GLuint attenuationTex;
 
+const float metaballScale = metaballR * 1.4f;
+
 void shadeClouds() {
   glDisable(GL_DEPTH_TEST);
   // shaderOutput * 0 + buffer * shader alpha
@@ -174,14 +180,17 @@ void shadeClouds() {
   glUniform1i(glGetUniformLocation(bbProg, "tex"), 0);
 
   GLuint modelLoc = glGetUniformLocation(bbProg, "model");
+  glUniform1i(glGetUniformLocation(bbProg, "debug"), 0);
 
   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);
+    mat4 model = translate(mat4(1), k.pos);
 
     // rotate the billboard so that its normal is oriented to the sun
     model = faceView(model);
 
+    model = scale(model, vec3(metaballScale));
+
     glUniformMatrix4fv(modelLoc, 1, GL_FALSE, glm::value_ptr(model));
 
     // Set the billboard color as RGBA = (1.0, 1.0, 1.0, 1.0).
@@ -213,7 +222,6 @@ void shadeClouds() {
                  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.
@@ -224,9 +232,20 @@ void shadeClouds() {
   checkError();
 }
 
-void renderObject() {}
+void renderObject() {
+  glDisable(GL_BLEND);
+  // render the sun
+  glUseProgram(sunProg);
+  mat4 model = translate(mat4(1), sunPos);
+  /* model = lookAt(sunPos, sunPos + sunDir, {0, 1, 0}) * model; */
+  model = translate(scale(translate(model, -sunPos), vec3(0.3)), sunPos);
+  glUniformMatrix4fv(glGetUniformLocation(sunProg, "model"), 1, GL_FALSE, glm::value_ptr(model));
+  glDrawElements(GL_TRIANGLES, 6, GL_UNSIGNED_INT, 0); 
+}
 
 void renderClouds() {
+  glUseProgram(bbProg);
+
   // Sort metaballs in descending order from the viewpoint
   sort(metaballs.begin(), metaballs.end(), [](Metaball &a, Metaball &b) {
     return distance(camPos, a.pos) > distance(camPos, b.pos);
@@ -248,10 +267,12 @@ void renderClouds() {
     GLuint modelLoc = glGetUniformLocation(bbProg, "model");
 
     // Place the billboard at the center of the corresponding metaball n.
-    mat4 model = scale(translate(mat4(1), k.pos), vec3(k.r) * 2.f);
+    mat4 model = translate(mat4(1), k.pos);
     // Rotate the billboard so that its normal is oriented to the viewpoint.
     model = faceView(model);
 
+    model = scale(model, vec3(metaballScale));
+
     glUniformMatrix4fv(modelLoc, 1, GL_FALSE, glm::value_ptr(model));
 
     // Set the billboard color as C[n].
@@ -266,6 +287,7 @@ void renderClouds() {
     // Don't modulate it -- blend it
     glUniform1i(glGetUniformLocation(bbProg, "modulate"), 0);
 
+    glUniform1f(glGetUniformLocation(bbProg, "debugColor"), curMode == debugColor);
     if (curMode != render) {
       float debugVal = 0;
       if (curMode == debugContDist) debugVal = k.d;
@@ -273,7 +295,7 @@ void renderClouds() {
       else if (curMode == debugProbExt) debugVal = cs.p_ext[k.coords.x][k.coords.y][k.coords.z];
       glUniform1f(glGetUniformLocation(bbProg, "debugVal"), debugVal);
       glDisable(GL_BLEND);
-      model = scale(model, vec3(0.02));
+      model = scale(model, vec3(0.2));
       glUniformMatrix4fv(modelLoc, 1, GL_FALSE, glm::value_ptr(model));
     }
 
@@ -298,6 +320,9 @@ void display() {
 
   view = glm::lookAt(camPos, viewPos, {0, 1, 0});
   proj = glm::perspective(45.f, aspect, znear, zfar);
+  glUseProgram(sunProg);
+  setProjectionAndViewUniforms(sunProg);
+  glUseProgram(bbProg);
   setProjectionAndViewUniforms(bbProg);
 
   glClearColor(0.83, 1, 1, 1); // background color
@@ -322,7 +347,7 @@ void keyboard(unsigned char key, int x, int y) {
   if (key == ' ') {
     calculateMetaballs();
     needsRedisplay = true;
-    needsReshading = true;
+    needsReshading = curMode == render;
   }
   if (key == '0') {
     curMode = render;
@@ -333,10 +358,14 @@ void keyboard(unsigned char key, int x, int y) {
     needsRedisplay = true;
   }
   if (key == '2') {
-    curMode = debugProbAct;
+    curMode = debugColor;
     needsRedisplay = true;
   }
   if (key == '3') {
+    curMode = debugProbAct;
+    needsRedisplay = true;
+  }
+  if (key == '4') {
     curMode = debugProbExt;
     needsRedisplay = true;
   }
@@ -353,7 +382,7 @@ void motion(int x, int y) {
   float dx = x - prevMouseX, dy = y - prevMouseY;
   prevMouseX = x;
   prevMouseY = y;
-  const vec3 origin(0, 18, 0);
+  const vec3 origin(0, 0, 0);
   const float sensitivity = 0.003f;
   auto camMat = translate(mat4(1), origin + camPos);
   auto rotation = rotate(rotate(mat4(1), -dx * sensitivity, {0, 1, 0}),
@@ -379,11 +408,14 @@ int main(int argc, char **argv) {
   glewInit();
 
   Program prog("billboardvert.glsl", "billboardfrag.glsl");
-
   bbProg = prog.progId;
-  glUseProgram(bbProg);
+  Program sProg("sunvert.glsl", "sunfrag.glsl");
+  sunProg = sProg.progId;
 
   glGenVertexArrays(1, &bbVao);
+  glUseProgram(sunProg);
+  glBindVertexArray(bbVao);
+  glUseProgram(bbProg);
   glBindVertexArray(bbVao);
   GLuint vbos[2];
   glGenBuffers(2, vbos);
index 360fb07481a2702188c6afd63ea9d34b72c27846..4e107e72670a2d86452d0383fcdf8563c786ce4a 100644 (file)
@@ -7,10 +7,10 @@ inline float randf() {
 }
 
 // Helper to account for bounds
-inline void set(float x[CLOUD_DIM][CLOUD_DIM][CLOUD_DIM], int i, int j, int k, float y) {
-  if (i < 0 || i >= CLOUD_DIM ||
-      j < 0 || j >= CLOUD_DIM ||
-      k < 0 || k >= CLOUD_DIM)
+inline void set(float x[CLOUD_DIM_X][CLOUD_DIM_Y][CLOUD_DIM_Z], int i, int j, int k, float y) {
+  if (i < 0 || i >= CLOUD_DIM_X ||
+      j < 0 || j >= CLOUD_DIM_Y ||
+      k < 0 || k >= CLOUD_DIM_Z)
     return;
   x[i][j][k] = y;
 }
@@ -20,12 +20,12 @@ inline void set(float x[CLOUD_DIM][CLOUD_DIM][CLOUD_DIM], int i, int j, int k, f
 #define P_ACT 0.001
 
 void initClouds(Clouds *cs) {
-  for (int i = 0; i < CLOUD_DIM; i++) {
-    for (int j = 0; j < CLOUD_DIM; j++) {
-      for (int k = 0; k < CLOUD_DIM; k++) {
+  for (int i = 0; i < CLOUD_DIM_X; i++) {
+    for (int j = 0; j < CLOUD_DIM_Y; j++) {
+      for (int k = 0; k < CLOUD_DIM_Z; k++) {
         cs->act[i][j][k] = randf() < 0.01;
         cs->cld[i][j][k] = false;
-        cs->hum[i][j][k] = randf() < 0.01;
+        cs->hum[i][j][k] = randf() < 0.1;
         cs->p_ext[i][j][k] = 0.f;
         cs->p_hum[i][j][k] = 0.f;
         cs->p_act[i][j][k] = 0.f;
@@ -33,22 +33,25 @@ void initClouds(Clouds *cs) {
     }
   }
 
-  for (int i = 0; i < CLOUD_DIM; i++) {
-    for (int j = 0; j < CLOUD_DIM; j++) {
-      for (int k = 0; k < CLOUD_DIM; k++) {
+  for (int i = 0; i < CLOUD_DIM_X; i++) {
+    for (int j = 0; j < CLOUD_DIM_Y; j++) {
+      for (int k = 0; k < CLOUD_DIM_Z; k++) {
         assert(cs->p_act[i][j][k] == 0.f);
         assert(cs->p_ext[i][j][k] == 0.f);
         assert(cs->p_hum[i][j][k] == 0.f);
       }
     }
   }
+  for (int k = 0; k < CLOUD_DIM_Z; k++)
+    cs->vz[k] = floor(randf() * 3);
 
   // generate ellipsoids of probability
-  for (int n = 0; n < 6; n++) {
+  const int numEllipsoids = CLOUD_DIM_X * CLOUD_DIM_Y * CLOUD_DIM_Z * 0.002;
+  for (int n = 0; n < numEllipsoids; n++) {
     const float maxSize = 8, minSize = 4;
     float delta = maxSize - minSize;
     int width = minSize + randf() * delta, height = minSize + randf() * delta, depth = minSize + randf() * delta;
-    int x = randf() * CLOUD_DIM, y = randf() * CLOUD_DIM, z = randf() * CLOUD_DIM;
+    int x = randf() * CLOUD_DIM_X, y = randf() * CLOUD_DIM_Y, z = randf() * CLOUD_DIM_Z;
     glm::vec3 center(x + width / 2, y + height / 2, z + depth / 2);
 
     for (int i = x; i < x + width; i++) {
@@ -65,10 +68,10 @@ void initClouds(Clouds *cs) {
 }
 
 // Helper to account for bounds
-inline bool get(bool x[CLOUD_DIM][CLOUD_DIM][CLOUD_DIM], int i, int j, int k) {
-  if (i < 0 || i >= CLOUD_DIM ||
-      j < 0 || j >= CLOUD_DIM ||
-      k < 0 || k >= CLOUD_DIM)
+inline bool get(bool x[CLOUD_DIM_X][CLOUD_DIM_Y][CLOUD_DIM_Z], int i, int j, int k) {
+  if (i < 0 || i >= CLOUD_DIM_X ||
+      j < 0 || j >= CLOUD_DIM_Y ||
+      k < 0 || k >= CLOUD_DIM_Z)
     return false;
   return x[i][j][k];
 }
@@ -83,9 +86,9 @@ inline bool f_act(Clouds *cs, int i, int j, int k) {
 void growth(Clouds *cs) {
   Clouds ncs = *cs;
 
-  for (int i = 0; i < CLOUD_DIM; i++) {
-    for (int j = 0; j < CLOUD_DIM; j++) {
-      for (int k = 0; k < CLOUD_DIM; k++) {
+  for (int i = 0; i < CLOUD_DIM_X; i++) {
+    for (int j = 0; j < CLOUD_DIM_Y; j++) {
+      for (int k = 0; k < CLOUD_DIM_Z; k++) {
         ncs.hum[i][j][k] = cs->hum[i][j][k] && !cs->act[i][j][k];
         ncs.cld[i][j][k] = cs->cld[i][j][k] || cs->act[i][j][k];
         ncs.act[i][j][k] = !cs->act[i][j][k] && cs->hum[i][j][k] && f_act(cs, i, j, k);
@@ -98,9 +101,9 @@ void growth(Clouds *cs) {
 
 void extinction(Clouds *cs) {
   Clouds ncs = *cs;
-  for (int i = 0; i < CLOUD_DIM; i++) {
-    for (int j = 0; j < CLOUD_DIM; j++) {
-      for (int k = 0; k < CLOUD_DIM; k++) {
+  for (int i = 0; i < CLOUD_DIM_X; i++) {
+    for (int j = 0; j < CLOUD_DIM_Y; j++) {
+      for (int k = 0; k < CLOUD_DIM_Z; k++) {
         ncs.cld[i][j][k] = cs->cld[i][j][k] && (randf() > cs->p_ext[i][j][k]);
         ncs.hum[i][j][k] = cs->hum[i][j][k] || (randf() < cs->p_hum[i][j][k]);
         ncs.act[i][j][k] = cs->act[i][j][k] || (randf() < cs->p_act[i][j][k]);
@@ -110,16 +113,34 @@ void extinction(Clouds *cs) {
   *cs = ncs;
 }
 
+void advection(Clouds *cs) {
+  Clouds ncs = *cs;
+
+  for (int i = 0; i < CLOUD_DIM_X; i++) {
+    for (int j = 0; j < CLOUD_DIM_Y; j++) {
+      for (int k = 0; k < CLOUD_DIM_Z; k++) {
+        int v = cs->vz[k];
+        ncs.hum[i][j][k] = i - v > 0 ? cs->hum[i - v][j][k] : 0;
+        ncs.cld[i][j][k] = i - v > 0 ? cs->cld[i - v][j][k] : 0;
+        ncs.act[i][j][k] = i - v > 0 ? cs->act[i - v][j][k] : 0;
+      }
+    }
+  }
+  
+  *cs = ncs;
+}
+
 /** Weighting function */
+// TODO: fill this out
 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++) {
+  for (int i = 0; i < CLOUD_DIM_X; i++) {
+    for (int j = 0; j < CLOUD_DIM_Y; j++) {
+      for (int k = 0; k < CLOUD_DIM_Z; k++) {
         float sum = 0;
 
         // sum
@@ -127,9 +148,9 @@ void calcContDist(Clouds *cls) {
           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)
+                if (i + ip < 0 || i + ip >= CLOUD_DIM_X ||
+                    j + jp < 0 || j + jp >= CLOUD_DIM_Y ||
+                    k + kp < 0 || k + kp >= CLOUD_DIM_Z)
                   continue;
 
                 sum += w(ip, jp, kp) * (float)cls->cld[i + ip][j + jp][k + kp];
@@ -148,5 +169,6 @@ void calcContDist(Clouds *cls) {
 void stepClouds(Clouds *cs) {
   growth(cs);
   extinction(cs);
+  advection(cs);
   calcContDist(cs);
 }
index c6c9e03b02cae0ae2938124685750bf7e8a793da..ef53425e5a8b99db86def28805b6ce886d9b3a93 100644 (file)
@@ -1,14 +1,21 @@
-#define CLOUD_DIM 16
+#define CLOUD_DIM_X 64
+#define CLOUD_DIM_Y 16
+#define CLOUD_DIM_Z 64
+/* #define CLOUD_DIM_X 16 */
+/* #define CLOUD_DIM_Y 16 */
+/* #define CLOUD_DIM_Z 16 */
 struct Clouds {
   // TODO: make more efficient
-  bool hum[CLOUD_DIM][CLOUD_DIM][CLOUD_DIM];
-  bool cld[CLOUD_DIM][CLOUD_DIM][CLOUD_DIM];
-  bool act[CLOUD_DIM][CLOUD_DIM][CLOUD_DIM];
-  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];
+  bool hum[CLOUD_DIM_X][CLOUD_DIM_Y][CLOUD_DIM_Z];
+  bool cld[CLOUD_DIM_X][CLOUD_DIM_Y][CLOUD_DIM_Z];
+  bool act[CLOUD_DIM_X][CLOUD_DIM_Y][CLOUD_DIM_Z];
+  float p_ext[CLOUD_DIM_X][CLOUD_DIM_Y][CLOUD_DIM_Z];
+  float p_hum[CLOUD_DIM_X][CLOUD_DIM_Y][CLOUD_DIM_Z];
+  float p_act[CLOUD_DIM_X][CLOUD_DIM_Y][CLOUD_DIM_Z];
+  /** Wind velocity */
+  float vz[CLOUD_DIM_Z];
   /** continuous distribution */
-  float q[CLOUD_DIM][CLOUD_DIM][CLOUD_DIM];
+  float q[CLOUD_DIM_X][CLOUD_DIM_Y][CLOUD_DIM_Z];
 };
 
 void initClouds(Clouds *cs);
diff --git a/sunfrag.glsl b/sunfrag.glsl
new file mode 100644 (file)
index 0000000..1d769fa
--- /dev/null
@@ -0,0 +1,5 @@
+#version 400
+out vec4 FragColor;
+void main() {
+  FragColor = vec4(1.f,1.f,1.f,1.f);
+}
diff --git a/sunvert.glsl b/sunvert.glsl
new file mode 100644 (file)
index 0000000..1aa37e1
--- /dev/null
@@ -0,0 +1,9 @@
+#version 400
+in vec3 vPosition;
+
+uniform mat4 model;
+uniform mat4 view;
+uniform mat4 projection;
+void main() {
+  gl_Position = projection * view * model * vec4(vPosition, 1);
+}