diff --git a/src/constants.h b/src/constants.h index 370d5b9..735f1b0 100755 --- a/src/constants.h +++ b/src/constants.h @@ -1,6 +1,9 @@ #pragma once #include "algebra.h" +#include +#include +#include /* ----------------------------------------------------------------------- | OPTIONS @@ -32,6 +35,12 @@ const static float DT = 0.0005f; // Time-step /* ----- GRID ----- */ const static float H_INV = 1.0f; // cell size ? +/* for sdf */ +const static glm::vec2 worldOrigin(0.f, 0.f); +const static float sdfCellSize = H_INV; +const static int sdfWidth = X_GRID; +const static int sdfHeight = Y_GRID; + /* ----- TRANSFER ----- */ #if INTERPOLATION == 1 const static int CUB = 2; diff --git a/src/main b/src/main index 5659430..3b62e7a 100755 Binary files a/src/main and b/src/main differ diff --git a/src/main.cpp b/src/main.cpp index c1cc25e..20d50c9 100755 --- a/src/main.cpp +++ b/src/main.cpp @@ -9,6 +9,7 @@ /* Declarations */ void initGLContext(); GLFWwindow *initGLFWContext(); +void drawSdf(GLFWwindow *); Solver *Simulation; int t_count = 0; int frameNumber = 0; @@ -36,9 +37,9 @@ void Initialization() { std::vector inParticles = Material::InitializeParticles(); std::vector inPolygons = Polygon::InitializePolygons(); - // std::cout << inPolygons.size() << '\n'; - Simulation = new Solver(inBorders, inNodes, inParticles, inPolygons); + + Simulation->computeSdf(); // for static objects, compute only once } void Update() { @@ -86,10 +87,11 @@ int main(int argc, char **argv) { #else /* [2] : show result on OpenGL window, and record an .mp4 if selected */ - const int frameLimit = 3; + const int frameLimit = 99999; GLFWwindow *window = initGLFWContext(); initGLContext(); + while (!glfwWindowShouldClose(window) && frameNumber < frameLimit) { glClear(GL_COLOR_BUFFER_BIT); @@ -99,42 +101,47 @@ int main(int argc, char **argv) { // AddParticles(); // P2G, compute grid forces, etc. - Update(); + // Update(); // Display frame at desired rate if (t_count % (int)(DT_render / DT) == 0) { Simulation->Draw(); - glfwSwapBuffers(window); - -#if RECORD_VIDEO - std::string dir = "../result/sim"; - // zero padding - std::string num = std::to_string(frameNumber); - num = std::string(4 - num.length(), '0') + num; - std::string output = dir + num + ".png"; + // for testing sdf + drawSdf(window); - // use opencv to save a frame - // the frame size becomes twice from frame 2, - // thus multiplying 2 is needed. - // don't know why - cv::Mat outImg(Y_WINDOW * 2, X_WINDOW * 2, CV_8UC3); - glPixelStorei(GL_PACK_ALIGNMENT, 4); - glReadPixels(0, 0, X_WINDOW * 2, Y_WINDOW * 2, GL_BGR, GL_UNSIGNED_BYTE, - (GLvoid *)outImg.data); - cv::flip(outImg, outImg, 0); - cv::imwrite(output, outImg); - - std::cout << "Frame " << frameNumber << " saved." << '\n'; + glfwSwapBuffers(window); - frameNumber++; -#endif + // #if RECORD_VIDEO + // std::string dir = "../result/sim"; + // + // // zero padding + // std::string num = std::to_string(frameNumber); + // num = std::string(4 - num.length(), '0') + num; + // std::string output = dir + num + ".png"; + // + // // use opencv to save a frame + // // the frame size becomes twice from frame 2, + // // thus multiplying 2 is needed. + // // don't know why + // cv::Mat outImg(Y_WINDOW * 2, X_WINDOW * 2, CV_8UC3); + // glPixelStorei(GL_PACK_ALIGNMENT, 4); + // glReadPixels(0, 0, X_WINDOW * 2, Y_WINDOW * 2, GL_BGR, + // GL_UNSIGNED_BYTE, + // (GLvoid *)outImg.data); + // cv::flip(outImg, outImg, 0); + // cv::imwrite(output, outImg); + // + // std::cout << "Frame " << frameNumber << " saved." << '\n'; + // + // frameNumber++; + // #endif glfwPollEvents(); } // end outer if // Don't forget to reset grid every frame !! - Simulation->ResetGrid(); + // Simulation->ResetGrid(); t_count++; } // end while @@ -187,3 +194,39 @@ void initGLContext() { glClearColor(.659f, .694f, .82f, .0f); glClear(GL_COLOR_BUFFER_BIT); } + +// for testing sdf +void drawSdf(GLFWwindow *wnd) { + double xpos, ypos; + glfwGetCursorPos(wnd, &xpos, &ypos); + + float fx, fy; + fx = xpos / (float)X_WINDOW; + fy = ypos / (float)Y_WINDOW; + + // Note: (0, 0) is at right-bottom + float wx, wy; + wx = fx * X_GRID; + wy = (1.f - fy) * Y_GRID; + + // for test + // draw sdf gradient at a given point + glm::vec2 start = glm::vec2(wx, wy); + float dist = Simulation->getDistance(start); + glm::vec2 grad = Simulation->getGradient(start); + glm::vec2 end = start + dist * grad; + + // start.x = start.x / (float)X_GRID * (float)X_WINDOW; + // start.y = start.y / (float)Y_GRID * (float)Y_WINDOW; + // + // end.x = end.x / (float)X_GRID * (float)X_WINDOW; + // end.y = end.y / (float)Y_GRID * (float)Y_WINDOW; + + glLineWidth(4); + glColor3f(1.f, 0.f, 0.f); + + glBegin(GL_LINES); + glVertex2f(start.x, start.y); // start point + glVertex2f(end.x, end.y); // end point + glEnd(); +} diff --git a/src/sdf.cpp b/src/sdf.cpp index df710ba..10dc048 100644 --- a/src/sdf.cpp +++ b/src/sdf.cpp @@ -1,8 +1,6 @@ #include "sdf.h" #include "solver.h" -glm::vec2 worldOrigin(0.f, 0.f); - /* Functions in class Polygon */ void Polygon::computeAabb() { lb = vertices[0]; @@ -97,116 +95,3 @@ void Polygon::DrawPolygon() { glEnd(); } /* end Functions in class Polygon */ - -bool intersect(glm::vec2 p1, glm::vec2 p2, glm::vec2 p3, glm::vec2 p4) { - glm::vec3 a(p1, 0), b(p2, 0), c(p3, 0), d(p4, 0); - - float crossAcAb, crossAbAd, crossDaDc, crossDcDb; - crossAcAb = glm::cross(c - a, b - a).z; - crossAbAd = glm::cross(b - a, d - a).z; - crossDaDc = glm::cross(a - d, c - d).z; - crossDcDb = glm::cross(c - d, b - d).z; - - if (crossAcAb * crossAbAd > 0 && crossDaDc * crossDcDb > 0) { - return true; - } else { - return false; - } -} - -bool inside_polygon(glm::vec2 p, Polygon &poly) { - int count = 0; - glm::vec2 q(1234567.f, 1234567.f); // a point at the infinity - int len = poly.vertices.size(); - - for (int i = 0; i < len; i++) { - glm::vec2 &start = poly.vertices[i]; - glm::vec2 &end = poly.vertices[(i + 1) % len]; - count += intersect(p, q, start, end); - } - - return count % 2 == 1; -} - -float nearest_distance(glm::vec2 p, Polygon &poly) { - float dist = 9999.f; - int len = poly.vertices.size(); - - for (int i = 0; i < len; i++) { - glm::vec2 &a = poly.vertices[i]; - glm::vec2 &b = poly.vertices[(i + 1) % len]; - - //当前处理的边定义为 ab - glm::vec2 ab = b - a; - float abLength = glm::length(ab); - glm::vec2 dir = glm::normalize(ab); - - //求点 p 在边 ab 上的投影点 q - //于是，线段 pq 的长度就是 p 点的有向距离 - glm::vec2 ap = p - a; - float frac = glm::dot(ap, dir); - frac = glm::clamp(frac, 0.f, abLength); - glm::vec2 q = a + frac * dir; - glm::vec2 pq = q - p; - - dist = glm::min(dist, length(pq)); - } - - return dist; -} - -// p is world position -// float getDistance(glm::vec2 p) { -// // world position to sdf index -// int idx_x, idx_y; -// idx_x = int(glm::floor(p.x / sdfCellSize)); -// idx_y = int(glm::floor(p.y / sdfCellSize)); -// -// if (idx_x < 0 || idx_x > sdfWidth - 1) { -// return 9999.f; -// } else if (idx_y < 0.f || idx_y > sdfHeight - 1) { -// return 9999.f; -// } else { -// return grid.sdf.at(Point(idx_x, idx_y)); -// } -// } - -// p is world position - -void computeSdf(int sdfWidth, int sdfHeight, float sdfCellSize, - std::vector &polygons, Solver &slv) { - // compute sdf for (sdfWidth * sdfHeight) world space points - // the interval between two adjacent points is sdfCellSize - for (int x = 0; x < sdfWidth; x++) { - for (int y = 0; y < sdfHeight; y++) { - // std::cout << x << ", " << y << '\n'; - - // world space point - glm::vec2 p = worldOrigin + glm::vec2(sdfCellSize * x, sdfCellSize * y); - - float fDist = 9999.f; - float temp = 0.f; - - for (int i = 0; i < polygons.size(); i++) { - // std::cout << i << '\n'; - temp = (inside_polygon(p, polygons[i])) ? -1.f : 1.f; - temp *= nearest_distance(p, polygons[i]); - // std::cout << "here" << '\n'; - - // fDist = glm::min(temp, fDist); - if (temp < fDist) { - fDist = temp; - - // change to Node - // the nearest polygon to this node - // grid.polyPtrs.set(y, x, polygons[i]); - } - } - - // save sdf distance - int idx = x + y * (sdfWidth - 1); // # of nodes in a row is (sdfWidth - 1) - slv.nodes[idx].sdfDist = fDist; - // std::cout << "(" << x << ", " << y << ") distance = " << fDist << '\n'; - } - } -} diff --git a/src/sdf.h b/src/sdf.h index 2320b1a..e07cf9e 100644 --- a/src/sdf.h +++ b/src/sdf.h @@ -10,8 +10,6 @@ #include #include -extern glm::vec2 worldOrigin; - class Polygon { public: std::vector vertices; @@ -51,7 +49,3 @@ class Polygon { return polys; } }; // end class Polygon - -bool intersect(glm::vec2, glm::vec2, glm::vec2, glm::vec2); -bool inside_polygon(glm::vec2, Polygon &); -float nearest_distance(glm::vec2, Polygon &); diff --git a/src/solver.cpp b/src/solver.cpp index d61a34c..5293763 100755 --- a/src/solver.cpp +++ b/src/solver.cpp @@ -16,6 +16,149 @@ Solver::Solver(const std::vector &inBorders, polylen = polygons.size(); } +/* Functions for computing sdf */ +void Solver::computeSdf() { + // compute sdf for (sdfWidth * sdfHeight) world space points + // the interval between two adjacent points is sdfCellSize + for (int x = 0; x < sdfWidth; x++) { + for (int y = 0; y < sdfHeight; y++) { + // std::cout << x << ", " << y << '\n'; + + // world space point + glm::vec2 p = worldOrigin + glm::vec2(sdfCellSize * x, sdfCellSize * y); + + float fDist = 9999.f; + float temp = 0.f; + + for (int i = 0; i < polygons.size(); i++) { + // std::cout << i << '\n'; + temp = (inside_polygon(p, polygons[i])) ? -1.f : 1.f; + temp *= nearest_distance(p, polygons[i]); + // std::cout << "here" << '\n'; + + // fDist = glm::min(temp, fDist); + if (temp < fDist) { + fDist = temp; + + // change to Node + // the nearest polygon to this node + // grid.polyPtrs.set(y, x, polygons[i]); + } + } + + // save sdf distance + int idx = x + y * (sdfWidth - 1); // # of nodes in a row is (sdfWidth - 1) + nodes[idx].sdfDist = fDist; + // std::cout << "(" << x << ", " << y << ") distance = " << fDist << '\n'; + } + } +} + +bool Solver::intersect(glm::vec2 p1, glm::vec2 p2, glm::vec2 p3, glm::vec2 p4) { + glm::vec3 a(p1, 0), b(p2, 0), c(p3, 0), d(p4, 0); + + float crossAcAb, crossAbAd, crossDaDc, crossDcDb; + crossAcAb = glm::cross(c - a, b - a).z; + crossAbAd = glm::cross(b - a, d - a).z; + crossDaDc = glm::cross(a - d, c - d).z; + crossDcDb = glm::cross(c - d, b - d).z; + + if (crossAcAb * crossAbAd > 0 && crossDaDc * crossDcDb > 0) { + return true; + } else { + return false; + } +} + +bool Solver::inside_polygon(glm::vec2 p, Polygon &poly) { + int count = 0; + glm::vec2 q(1234567.f, 1234567.f); // a point at the infinity + int len = poly.vertices.size(); + + for (int i = 0; i < len; i++) { + glm::vec2 &start = poly.vertices[i]; + glm::vec2 &end = poly.vertices[(i + 1) % len]; + count += intersect(p, q, start, end); + } + + return count % 2 == 1; +} + +float Solver::nearest_distance(glm::vec2 p, Polygon &poly) { + float dist = 9999.f; + int len = poly.vertices.size(); + + for (int i = 0; i < len; i++) { + glm::vec2 &a = poly.vertices[i]; + glm::vec2 &b = poly.vertices[(i + 1) % len]; + + //当前处理的边定义为 ab + glm::vec2 ab = b - a; + float abLength = glm::length(ab); + glm::vec2 dir = glm::normalize(ab); + + //求点 p 在边 ab 上的投影点 q + //于是，线段 pq 的长度就是 p 点的有向距离 + glm::vec2 ap = p - a; + float frac = glm::dot(ap, dir); + frac = glm::clamp(frac, 0.f, abLength); + glm::vec2 q = a + frac * dir; + glm::vec2 pq = q - p; + + dist = glm::min(dist, length(pq)); + } + + return dist; +} + +// p is world position +float Solver::getDistance(glm::vec2 p) { + + // world position to sdf index + int idx_x, idx_y; + idx_x = int(glm::floor(p.x / sdfCellSize)); + idx_y = int(glm::floor(p.y / sdfCellSize)); + + if (idx_x < 0 || idx_x > sdfWidth - 1) { + return 9999.f; + } else if (idx_y < 0.f || idx_y > sdfHeight - 1) { + return 9999.f; + } else { + int idx = idx_x + idx_y * (sdfWidth - 1); + return nodes[idx].sdfDist; + } +} + +// p is world position +glm::vec2 Solver::getGradient(glm::vec2 p) { + glm::vec2 gd; + + // world position to sdf position + glm::vec2 sdfPos, sdfPosFloor; + sdfPos = p / sdfCellSize; + sdfPosFloor = glm::floor(sdfPos); + + float fx, fy; + fx = sdfPos.x - sdfPosFloor.x; + fy = sdfPos.y - sdfPosFloor.y; + + float temp1 = getDistance(glm::vec2(p.x + 1.f, p.y)) - + getDistance(glm::vec2(p.x - 1.f, p.y)); + float temp2 = getDistance(glm::vec2(p.x + 1.f, p.y + 1.f)) - + getDistance(glm::vec2(p.x - 1.f, p.y + 1.f)); + gd.x = glm::lerp(temp1, temp2, fy); + + float temp3 = getDistance(glm::vec2(p.x, p.y + 1.f)) - + getDistance(glm::vec2(p.x, p.y - 1.f)); + float temp4 = getDistance(glm::vec2(p.x + 1.f, p.y + 1.f)) - + getDistance(glm::vec2(p.x + 1.f, p.y - 1.f)); + gd.y = glm::lerp(temp3, temp4, fx); + + gd = glm::normalize(-gd); + + return gd; +} + /* ----------------------------------------------------------------------- | MATERIAL POINT METHOD ALGORITHM | diff --git a/src/solver.h b/src/solver.h index 264de7f..9a2b86d 100755 --- a/src/solver.h +++ b/src/solver.h @@ -2,6 +2,7 @@ #include "node.h" #include "particle.h" +#include "constants.h" #include "sdf.h" #include #include @@ -35,6 +36,14 @@ class Solver { void UpdateParticles(); void ResetGrid(); + /* Functions for computing sdf */ + void computeSdf(); + bool intersect(glm::vec2, glm::vec2, glm::vec2, glm::vec2); + bool inside_polygon(glm::vec2, Polygon &); + float nearest_distance(glm::vec2, Polygon &); + float getDistance(glm::vec2); + glm::vec2 getGradient(glm::vec2); + void Draw(); // Draw particles, border and nodes (if selected) void WriteToFile( int frame); // Write point cloud coordinates to .ply file (Houdini)