From 0258cd283dc666dd89f620179f57343662ab5d1a Mon Sep 17 00:00:00 2001 From: Luke Lau Date: Fri, 10 Apr 2020 15:02:27 +0100 Subject: [PATCH] Continuous distribution --- .gitignore | 1 + Makefile | 2 +- billboardfrag.glsl | 6 ++ clouds.cpp | 194 +++++++++++++++++++++++---------------------- debug.hpp | 31 +++++++- simulation.cpp | 36 +++++++++ simulation.h | 3 +- 7 files changed, 173 insertions(+), 100 deletions(-) diff --git a/.gitignore b/.gitignore index dfe6a70..309bd5c 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ clouds fbo.tga *.dSYM +bbtex diff --git a/Makefile b/Makefile index 67b688d..6465814 100644 --- 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 diff --git a/billboardfrag.glsl b/billboardfrag.glsl index 94d15dc..046b7bd 100644 --- a/billboardfrag.glsl +++ b/billboardfrag.glsl @@ -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 diff --git a/clouds.cpp b/clouds.cpp index d5c8f6a..52b3053 100644 --- a/clouds.cpp +++ b/clouds.cpp @@ -1,27 +1,26 @@ #include "debug.hpp" -#include "simulation.h" #include "program.hpp" +#include "simulation.h" #include #include +#include #include #include #include #include #include +#include #include #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 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 metaballs = {{{-0.5, 0.5, 0.5}, 0.25}, */ -/* {{-0.3, 0.5, 0.3}, 0.25}}; */ -vector metaballs = {{{0, 0, 0.5}, 1.f}, - {{0, 0.3, 0.3}, 0.7f}}; + +array 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 diff --git a/debug.hpp b/debug.hpp index debe65e..6774a1a 100644 --- 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); +} diff --git a/simulation.cpp b/simulation.cpp index 0471f4a..e022d12 100644 --- a/simulation.cpp +++ b/simulation.cpp @@ -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); } diff --git a/simulation.h b/simulation.h index f021ce4..c6c9e03 100644 --- a/simulation.h +++ b/simulation.h @@ -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); -- 2.30.2