From ce1a6f5e411a92726fc1887f604fc18df15ebe01 Mon Sep 17 00:00:00 2001 From: Jiaming Guo Date: Fri, 11 Aug 2017 09:13:17 -0500 Subject: [PATCH] add functions for incorporating nwt files and prepare for computing largest connected component --- flow.h | 537 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- main.cu | 165 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------- 2 files changed, 416 insertions(+), 286 deletions(-) diff --git a/flow.h b/flow.h index b1fa43d..dc7890a 100644 --- a/flow.h +++ b/flow.h @@ -355,6 +355,24 @@ namespace stim { return E[i].r(j); } + // get the pendant vertices + std::vector get_pendant_vertex() { + std::vector result; + int count = 0; + + for (unsigned i = 0; i < V.size(); i++) { // for every vertex + for (unsigned j = 0; j < E.size(); j++) { // for every edge + if (i == E[j].v[0] || i == E[j].v[1]) // check whether current vertex terminates one edge + count++; + } + if (count == 1) // is pendant vertex + result.push_back(i); + count = 0; // reset count + } + + return result; + } + // get the velocity of pendant vertex i T get_velocity(unsigned i) { @@ -378,6 +396,12 @@ namespace stim { P[i] = value; } + // extrct the largest connected component + void extract_lcc() { + + + } + // solve the linear system to get stable flow state void solve_flow(T viscosity) { @@ -385,7 +409,7 @@ namespace stim { clear(); // get the pendant vertex indices - pendant_vertex = get_boundary_vertex(); + pendant_vertex = get_pendant_vertex(); // get bounding box bb = (*this).boundingbox(); @@ -749,7 +773,7 @@ namespace stim { } // draw solid sphere at every vertex - void glSolidSphere(T max_pressure, T scale, GLint subdivision) { + void glSolidSphere(T max_pressure, GLint subdivision, T scale = 1.0f) { // waste? for (unsigned i = 0; i < num_edge; i++) { @@ -808,7 +832,7 @@ namespace stim { } // draw edges as series of cylinders - void glSolidCylinder(unsigned index, std::vector color, T scale, GLint subdivision) { + void glSolidCylinder(unsigned index, std::vector color, GLint subdivision, T scale = 1.0f) { stim::vec3 tmp_d; stim::vec3 center1; @@ -895,7 +919,7 @@ namespace stim { } // draw the flow direction as cone, the size of the cone depends on the length of that edge - void glSolidCone(unsigned i, T scale, GLint subdivision) { + void glSolidCone(unsigned i, GLint subdivision, T scale = 1.0f, T threshold = 0.01f) { stim::vec3 tmp_d; // direction stim::vec3 center; // cone hat center @@ -909,27 +933,28 @@ namespace stim { unsigned index = E[i].size() / 2 - 1; tmp_d = E[i][index + 1] - E[i][index]; - h = tmp_d.len() / 3.0f; // get the height base by factor 3 + h = tmp_d.len() / 1.5f; // get the height base by factor 3 tmp_d = tmp_d.norm(); center = (E[i][index + 1] + E[i][index]) / 2; tmp_c.rotate(tmp_d); radius = (E[i].r(index + 1) + E[i].r(index)) / 2 * scale; radius = (h / sqrt(3) < radius) ? h / sqrt(3) : radius; // update radius - if (v[i] > 0) + if (v[i] > threshold) head = center + tmp_d * h; - else + else if (v[i] < -threshold) head = center - tmp_d * h; stim::circle c(center, radius, tmp_d, tmp_c.U); cp = c.glpoints(subdivision); - glBegin(GL_TRIANGLE_FAN); - glVertex3f(head[0], head[1], head[2]); - for (unsigned k = 0; k < cp.size(); k++) - glVertex3f(cp[k][0], cp[k][1], cp[k][2]); - glEnd(); - glFlush(); - + if (v[i] > threshold || v[i] < -threshold) { + glBegin(GL_TRIANGLE_FAN); + glVertex3f(head[0], head[1], head[2]); + for (unsigned k = 0; k < cp.size(); k++) + glVertex3f(cp[k][0], cp[k][1], cp[k][2]); + glEnd(); + glFlush(); + } // draw a cone for every edge to indicate //for (unsigned j = 0; j < E[i].size() - 1; j++) { // for every point on current edge // tmp_d = E[i][j + 1] - E[i][j]; @@ -953,7 +978,7 @@ namespace stim { //} //glFlush(); } - void glSolidCone(GLint subdivision) { + void glSolidCone(GLint subdivision, T scale = 1.0f, T threhold = 0.01f) { stim::vec3 tmp_d; // direction stim::vec3 center; // cone hat center @@ -973,21 +998,22 @@ namespace stim { tmp_d = tmp_d.norm(); center = (E[i][k2] + E[i][k1]) / 2; tmp_c.rotate(tmp_d); - radius = (E[i].r(k2) + E[i].r(k1)) / 2; + radius = (E[i].r(k2) + E[i].r(k1)) / 2 * scale; radius = (h / sqrt(3) < radius) ? h / sqrt(3) : radius; // update radius by height base if necessary - if (v[i] > 0) // if flow flows from k1 to k2 + if (v[i] > threhold) // if flow flows from k1 to k2 head = center + tmp_d * h; - else + else if(v[i] < -threhold) head = center - tmp_d * h; stim::circle c(center, radius, tmp_d, tmp_c.U); cp = c.glpoints(subdivision); - glBegin(GL_TRIANGLE_FAN); - glVertex3f(head[0], head[1], head[2]); - for (unsigned k = 0; k < cp.size(); k++) - glVertex3f(cp[k][0], cp[k][1], cp[k][2]); - glEnd(); - + if (v[i] > threhold || v[i] < -threhold) { + glBegin(GL_TRIANGLE_FAN); + glVertex3f(head[0], head[1], head[2]); + for (unsigned k = 0; k < cp.size(); k++) + glVertex3f(cp[k][0], cp[k][1], cp[k][2]); + glEnd(); + } //for (unsigned j = 0; j < E[i].size() - 1; j++) { // for every point on current edge // tmp_d = E[i][j + 1] - E[i][j]; // tmp_d = tmp_d.norm(); @@ -1013,7 +1039,7 @@ namespace stim { } // draw main feeder as solid cube - void glSolidCuboid(bool arrow, GLint subdivision, bool manufacture = false, T length = 40.0f, T height = 10.0f) { + void glSolidCuboid(GLint subdivision, bool manufacture = false, T length = 40.0f, T height = 10.0f) { T width; stim::vec3 L = bb.A; // get the bottom left corner @@ -1073,51 +1099,10 @@ namespace stim { glVertex3f(main_feeder[i][0] + length / 2, main_feeder[i][1] - height / 2, main_feeder[i][2] + width / 2); glEnd(); } - - // render total flow direction - if (arrow) { - stim::vec3 d = main_feeder[1] - main_feeder[0]; - d = d.norm(); - stim::vec3 v1, v2, v3; - stim::vec3 center = bb.center(); - stim::vec3 size = bb.size(); - stim::circle tmp_c; - tmp_c.rotate(d); - std::vector > cp1(subdivision + 1); - std::vector > cp2(subdivision + 1); - v2 = center - stim::vec3(0.0f, size[1] + 10.0f, 0.0f); - v1 = v2 - stim::vec3(30.0f, 0.0f, 0.0f); - v3 = v2 + stim::vec3(30.0f, 0.0f, 0.0f); - v2 = v2 + stim::vec3(20.0f, 0.0f, 0.0f); - - // draw the total flow direciton indicating arrow - stim::circle c1(v1, 5.0f / 2, d, tmp_c.U); - cp1 = c1.glpoints(subdivision); - stim::circle c2(v2, 5.0f / 2, d, tmp_c.U); - cp2 = c2.glpoints(subdivision); - - glColor3f(0.0f, 1.0f, 0.0f); - glBegin(GL_QUAD_STRIP); - for (unsigned k = 0; k < cp1.size(); k++) { - glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]); - glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]); - } - glEnd(); - - // render the cone part - stim::circle c3(v2, 5.0f, d, tmp_c.U); - cp2 = c3.glpoints(subdivision); - glBegin(GL_TRIANGLE_FAN); - glVertex3f(v3[0], v3[1], v3[2]); - for (unsigned k = 0; k < cp2.size(); k++) - glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]); - glEnd(); - } - glFlush(); } // draw flow velocity field, glyph - void glyph(std::vector color, GLint subdivision, T r = 4.0f) { + void glyph(std::vector color, GLint subdivision, T scale = 1.0f, bool frame = false, T r = 4.0f, T threshold = 0.01f) { // v1----v2-->v3 T k = 4.0f; // quartering @@ -1127,51 +1112,124 @@ namespace stim { std::vector > cp1(subdivision + 1); std::vector > cp2(subdivision + 1); + // rendering the arrows for (unsigned i = 0; i < num_edge; i++) { // for every edge glColor3f((float)color[i * 3 + 0] / 255.0f, (float)color[i * 3 + 1] / 255.0f, (float)color[i * 3 + 2] / 255.0f); for (unsigned j = 0; j < E[i].size() - 1; j++) { // for every point on that edge // consider the velocity valuence - if (v[i] > 0) { // positive, from start point to end point + if (v[i] > threshold) { // positive, from start point to end point v1 = E[i][j]; v3 = E[i][j + 1]; } - else { // negative, from end point to start point + else if (v[i] < -threshold) { // negative, from end point to start point v1 = E[i][j + 1]; v3 = E[i][j]; } - d = v3 - v1; - // place the arrow in the middel of one edge - v2 = v1 + (1.0f / k * 2.0f) * d; // looks like =->= - v1 = v1 + (1.0f / k) * d; - v3 = v3 - (1.0f / k) * d; - d = d.norm(); - tmp_c.rotate(d); - - // render the cylinder part - stim::circle c1(v1, r / 2, d, tmp_c.U); - cp1 = c1.glpoints(subdivision); - stim::circle c2(v2, r / 2, d, tmp_c.U); - cp2 = c2.glpoints(subdivision); - - glBegin(GL_QUAD_STRIP); - for (unsigned k = 0; k < cp1.size(); k++) { - glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]); - glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]); + + if (v[i] > threshold || v[i] < -threshold) { + d = v3 - v1; + // place the arrow in the middel of one edge + v2 = v1 + (1.0f / k * 2.0f) * d; // looks like =->= + v1 = v1 + (1.0f / k) * d; + v3 = v3 - (1.0f / k) * d; + d = d.norm(); + tmp_c.rotate(d); + + // render the cylinder part + stim::circle c1(v1, r / 2 * scale, d, tmp_c.U); + cp1 = c1.glpoints(subdivision); + stim::circle c2(v2, r / 2 * scale, d, tmp_c.U); + cp2 = c2.glpoints(subdivision); + + glBegin(GL_QUAD_STRIP); + for (unsigned k = 0; k < cp1.size(); k++) { + glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]); + glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]); + } + glEnd(); + + // render the cone part + stim::circle c3(v2, r * scale, d, tmp_c.U); + cp2 = c3.glpoints(subdivision); + glBegin(GL_TRIANGLE_FAN); + glVertex3f(v3[0], v3[1], v3[2]); + for (unsigned k = 0; k < cp2.size(); k++) + glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]); + glEnd(); } - glEnd(); + } + } - // render the cone part - stim::circle c3(v2, r, d, tmp_c.U); - cp2 = c3.glpoints(subdivision); - glBegin(GL_TRIANGLE_FAN); - glVertex3f(v3[0], v3[1], v3[2]); - for (unsigned k = 0; k < cp2.size(); k++) - glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]); - glEnd(); + // rendering frames + if (frame) { + frame = false; + stim::vec3 center1; + stim::vec3 center2; + stim::vec3 tmp_d; // flow direction + float r1, r2; + + for (unsigned i = 0; i < num_edge; i++) { + for (unsigned j = 0; j < E[i].size() - 1; j++) { + center1 = E[i][j]; + center2 = E[i][j + 1]; + + r1 = get_r(i, j) * scale; + r2 = get_r(i, j + 1) * scale; + + subdivision = 5; // rough frames + + if (j == 0) { + if (E[i].size() == 2) + find_envelope(cp1, cp2, center1, center2, r1, r2, subdivision); + else { + tmp_d = center2 - center1; + tmp_d = tmp_d.norm(); + tmp_c.rotate(tmp_d); + stim::circle c1(center1, r1, tmp_d, tmp_c.U); + cp1 = c1.glpoints(subdivision); + tmp_d = (E[i][j + 2] - center2) + (center2 - center1); + tmp_d = tmp_d.norm(); + tmp_c.rotate(tmp_d); + stim::circle c2(center2, r2, tmp_d, tmp_c.U); + cp2 = c2.glpoints(subdivision); + } + } + else if (j == E[i].size() - 2) { + tmp_d = (center2 - center1) + (center1 - E[i][j - 1]); + tmp_d = tmp_d.norm(); + tmp_c.rotate(tmp_d); + stim::circle c1(center1, r1, tmp_d, tmp_c.U); + cp1 = c1.glpoints(subdivision); + tmp_d = center2 - center1; + tmp_d = tmp_d.norm(); + tmp_c.rotate(tmp_d); + stim::circle c2(center2, r2, tmp_d, tmp_c.U); + cp2 = c2.glpoints(subdivision); + } + else { + tmp_d = (center2 - center1) + (center1 - E[i][j - 1]); + tmp_d = tmp_d.norm(); + tmp_c.rotate(tmp_d); + stim::circle c1(center1, r1, tmp_d, tmp_c.U); + cp1 = c1.glpoints(subdivision); + tmp_d = (E[i][j + 2] - center2) + (center2 - center1); + tmp_d = tmp_d.norm(); + tmp_c.rotate(tmp_d); + stim::circle c2(center2, r2, tmp_d, tmp_c.U); + cp2 = c2.glpoints(subdivision); + } + + glColor3f(140/255.0f, 81/255.0f, 10/255.0f); + glBegin(GL_LINES); + for (unsigned k = 0; k < cp1.size(); k++) { + glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]); + glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]); + } + glEnd(); + } } } - glFlush(); } // display the total volume flow rate @@ -1204,145 +1262,121 @@ namespace stim { } // draw the bridge as lines or arrows - void line_bridge(bool &redisplay, bool arrow, T r = 4.0f) { + void line_bridge(bool &redisplay, T r = 4.0f) { - if (redisplay) + if (redisplay) { // check to see whether the display list needs to be updated glDeleteLists(dlist, 1); - redisplay = false; + redisplay = false; + } if (!glIsList(dlist)) { dlist = glGenLists(1); glNewList(dlist, GL_COMPILE); - // render flow direction arrows - if (arrow) { - // v1----v2-->v3 - T k = 4.0f; // quartering - stim::vec3 v1, v2, v3; // three point - stim::vec3 d; // direction vector - stim::circle tmp_c; - std::vector > cp1(subdivision + 1); - std::vector > cp2(subdivision + 1); - - // inlet, right-going - for (unsigned i = 0; i < inlet.size(); i++) { - if (inlet_feasibility[i]) - glColor3f(0.0f, 0.0f, 0.0f); // feasible -> black - else - glColor3f(1.0f, 0.0f, 0.0f); // nonfeasible -> red - for (unsigned j = 0; j < inlet[i].V.size() - 1; j++) { - v1 = inlet[i].V[j]; - v3 = inlet[i].V[j + 1]; - d = v3 - v1; - // place the arrow in the middel of one edge - v2 = v1 + (1.0f / k * 2.0f) * d; // looks like =->= - v1 = v1 + (1.0f / k) * d; - v3 = v3 - (1.0f / k) * d; - d = d.norm(); - tmp_c.rotate(d); - - // render the cylinder part - stim::circle c1(v1, r / 2, d, tmp_c.U); - cp1 = c1.glpoints(subdivision); - stim::circle c2(v2, r / 2, d, tmp_c.U); - cp2 = c2.glpoints(subdivision); - - glBegin(GL_QUAD_STRIP); - for (unsigned k = 0; k < cp1.size(); k++) { - glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]); - glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]); - } - glEnd(); - - // render the cone part - stim::circle c3(v2, r, d, tmp_c.U); - cp2 = c3.glpoints(subdivision); - glBegin(GL_TRIANGLE_FAN); - glVertex3f(v3[0], v3[1], v3[2]); - for (unsigned k = 0; k < cp2.size(); k++) - glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]); - glEnd(); - } - } - - // outlet, right-going - for (unsigned i = 0; i < outlet.size(); i++) { - if (outlet_feasibility[i]) - glColor3f(0.0f, 0.0f, 0.0f); // feasible -> black - else - glColor3f(1.0f, 0.0f, 0.0f); // nonfeasible -> red - for (unsigned j = 0; j < outlet[i].V.size() - 1; j++) { - v1 = outlet[i].V[j + 1]; - v3 = outlet[i].V[j]; - d = v3 - v1; - // place the arrow in the middel of one edge - v2 = v1 + (1.0f / k * 2.0f) * d; // looks like =->= - v1 = v1 + (1.0f / k) * d; - v3 = v3 - (1.0f / k) * d; - d = d.norm(); - tmp_c.rotate(d); - - // render the cylinder part - stim::circle c1(v1, r / 2, d, tmp_c.U); - cp1 = c1.glpoints(subdivision); - stim::circle c2(v2, r / 2, d, tmp_c.U); - cp2 = c2.glpoints(subdivision); - - glBegin(GL_QUAD_STRIP); - for (unsigned k = 0; k < cp1.size(); k++) { - glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]); - glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]); - } - glEnd(); - - // render the cone part - stim::circle c3(v2, r, d, tmp_c.U); - cp2 = c3.glpoints(subdivision); - glBegin(GL_TRIANGLE_FAN); - glVertex3f(v3[0], v3[1], v3[2]); - for (unsigned k = 0; k < cp2.size(); k++) - glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]); - glEnd(); - } - } - - // render transparent lines as indexing - glEnable(GL_BLEND); // enable color blend - glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); // set blend function - glDisable(GL_DEPTH_TEST); - glLineWidth(5); - for (unsigned i = 0; i < inlet.size(); i++) { - if (inlet_feasibility[i]) - glColor4f(0.0f, 0.0f, 0.0f, 0.2f); // feasible -> black - else - glColor4f(1.0f, 0.0f, 0.0f, 0.2f); // nonfeasible -> red - - glBegin(GL_LINE_STRIP); - for (unsigned j = 0; j < inlet[i].V.size(); j++) - glVertex3f(inlet[i].V[j][0], inlet[i].V[j][1], inlet[i].V[j][2]); - glEnd(); - } - for (unsigned i = 0; i < outlet.size(); i++) { - if (outlet_feasibility[i]) - glColor4f(0.0f, 0.0f, 0.0f, 0.2f); // feasible -> black - else - glColor4f(1.0f, 0.0f, 0.0f, 0.2f); // nonfeasible -> red - glBegin(GL_LINE_STRIP); - for (unsigned j = 0; j < outlet[i].V.size(); j++) - glVertex3f(outlet[i].V[j][0], outlet[i].V[j][1], outlet[i].V[j][2]); - glEnd(); - } - glDisable(GL_BLEND); - glEnable(GL_DEPTH_TEST); - } - // render connection lines - else { + //// render flow direction arrows + //if (arrow) { + // // v1----v2-->v3 + // T k = 4.0f; // quartering + // stim::vec3 v1, v2, v3; // three point + // stim::vec3 d; // direction vector + // stim::circle tmp_c; + // std::vector > cp1(subdivision + 1); + // std::vector > cp2(subdivision + 1); + + // // inlet, right-going + // for (unsigned i = 0; i < inlet.size(); i++) { + // if (inlet_feasibility[i]) + // glColor3f(0.0f, 0.0f, 0.0f); // feasible -> black + // else + // glColor3f(1.0f, 0.0f, 0.0f); // nonfeasible -> red + // for (unsigned j = 0; j < inlet[i].V.size() - 1; j++) { + // v1 = inlet[i].V[j]; + // v3 = inlet[i].V[j + 1]; + // d = v3 - v1; + // // place the arrow in the middel of one edge + // v2 = v1 + (1.0f / k * 2.0f) * d; // looks like =->= + // v1 = v1 + (1.0f / k) * d; + // v3 = v3 - (1.0f / k) * d; + // d = d.norm(); + // tmp_c.rotate(d); + + // // render the cylinder part + // stim::circle c1(v1, r / 2, d, tmp_c.U); + // cp1 = c1.glpoints(subdivision); + // stim::circle c2(v2, r / 2, d, tmp_c.U); + // cp2 = c2.glpoints(subdivision); + + // glBegin(GL_QUAD_STRIP); + // for (unsigned k = 0; k < cp1.size(); k++) { + // glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]); + // glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]); + // } + // glEnd(); + + // // render the cone part + // stim::circle c3(v2, r, d, tmp_c.U); + // cp2 = c3.glpoints(subdivision); + // glBegin(GL_TRIANGLE_FAN); + // glVertex3f(v3[0], v3[1], v3[2]); + // for (unsigned k = 0; k < cp2.size(); k++) + // glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]); + // glEnd(); + // } + // } + + // // outlet, right-going + // for (unsigned i = 0; i < outlet.size(); i++) { + // if (outlet_feasibility[i]) + // glColor3f(0.0f, 0.0f, 0.0f); // feasible -> black + // else + // glColor3f(1.0f, 0.0f, 0.0f); // nonfeasible -> red + // for (unsigned j = 0; j < outlet[i].V.size() - 1; j++) { + // v1 = outlet[i].V[j + 1]; + // v3 = outlet[i].V[j]; + // d = v3 - v1; + // // place the arrow in the middel of one edge + // v2 = v1 + (1.0f / k * 2.0f) * d; // looks like =->= + // v1 = v1 + (1.0f / k) * d; + // v3 = v3 - (1.0f / k) * d; + // d = d.norm(); + // tmp_c.rotate(d); + + // // render the cylinder part + // stim::circle c1(v1, r / 2, d, tmp_c.U); + // cp1 = c1.glpoints(subdivision); + // stim::circle c2(v2, r / 2, d, tmp_c.U); + // cp2 = c2.glpoints(subdivision); + + // glBegin(GL_QUAD_STRIP); + // for (unsigned k = 0; k < cp1.size(); k++) { + // glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]); + // glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]); + // } + // glEnd(); + + // // render the cone part + // stim::circle c3(v2, r, d, tmp_c.U); + // cp2 = c3.glpoints(subdivision); + // glBegin(GL_TRIANGLE_FAN); + // glVertex3f(v3[0], v3[1], v3[2]); + // for (unsigned k = 0; k < cp2.size(); k++) + // glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]); + // glEnd(); + // } + // } + + // // render transparent lines as indexing + // glEnable(GL_BLEND); // enable color blend + // glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); // set blend function + // glDisable(GL_DEPTH_TEST); glLineWidth(5); for (unsigned i = 0; i < inlet.size(); i++) { if (inlet_feasibility[i]) - glColor3f(0.0f, 0.0f, 0.0f); // feasible -> black + glColor3f(0.0f, 0.0f, 0.0f); + // glColor4f(0.0f, 0.0f, 0.0f, 0.2f); // feasible -> black else - glColor3f(1.0f, 0.0f, 0.0f); // nonfeasible -> red + glColor3f(1.0f, 0.0f, 0.0f); + // glColor4f(1.0f, 0.0f, 0.0f, 0.2f); // nonfeasible -> red glBegin(GL_LINE_STRIP); for (unsigned j = 0; j < inlet[i].V.size(); j++) @@ -1351,23 +1385,31 @@ namespace stim { } for (unsigned i = 0; i < outlet.size(); i++) { if (outlet_feasibility[i]) - glColor3f(0.0f, 0.0f, 0.0f); // feasible -> black + glColor3f(0.0f, 0.0f, 0.0f); + // glColor4f(0.0f, 0.0f, 0.0f, 0.2f); // feasible -> black else - glColor3f(1.0f, 0.0f, 0.0f); // nonfeasible -> red + glColor3f(1.0f, 0.0f, 0.0f); + // glColor4f(1.0f, 0.0f, 0.0f, 0.2f); // nonfeasible -> red glBegin(GL_LINE_STRIP); for (unsigned j = 0; j < outlet[i].V.size(); j++) glVertex3f(outlet[i].V[j][0], outlet[i].V[j][1], outlet[i].V[j][2]); glEnd(); } - } - glFlush(); + // glDisable(GL_BLEND); + // glEnable(GL_DEPTH_TEST); + //} glEndList(); } glCallList(dlist); } // draw the bridge as tubes - void tube_bridge(T subdivision, T radius = 5.0f) { + void tube_bridge(bool &redisplay, T subdivision, T scale = 1.0f, T radius = 5.0f) { + + if (redisplay) { + glDeleteLists(dlist + 1, 1); + redisplay = false; + } if (!glIsList(dlist + 1)) { glNewList(dlist + 1, GL_COMPILE); @@ -1383,7 +1425,7 @@ namespace stim { for (unsigned j = 1; j < inlet[i].V.size() - 1; j++) { glPushMatrix(); glTranslatef(inlet[i].V[j][0], inlet[i].V[j][1], inlet[i].V[j][2]); - glutSolidSphere(radius, subdivision, subdivision); + glutSolidSphere(radius * scale, subdivision, subdivision); glPopMatrix(); } // render edge as cylinder @@ -1391,8 +1433,8 @@ namespace stim { dir = inlet[i].V[j] - inlet[i].V[j + 1]; dir = dir.norm(); unit_c.rotate(dir); - stim::circle c1(inlet[i].V[j], inlet[i].r, dir, unit_c.U); - stim::circle c2(inlet[i].V[j + 1], inlet[i].r, dir, unit_c.U); + stim::circle c1(inlet[i].V[j], inlet[i].r * scale, dir, unit_c.U); + stim::circle c2(inlet[i].V[j + 1], inlet[i].r * scale, dir, unit_c.U); cp1 = c1.glpoints(subdivision); cp2 = c2.glpoints(subdivision); @@ -1410,7 +1452,7 @@ namespace stim { for (unsigned j = 1; j < outlet[i].V.size() - 1; j++) { glPushMatrix(); glTranslatef(outlet[i].V[j][0], outlet[i].V[j][1], outlet[i].V[j][2]); - glutSolidSphere(radius, subdivision, subdivision); + glutSolidSphere(radius * scale, subdivision, subdivision); glPopMatrix(); } // render edge as cylinder @@ -1418,8 +1460,8 @@ namespace stim { dir = outlet[i].V[j] - outlet[i].V[j + 1]; dir = dir.norm(); unit_c.rotate(dir); - stim::circle c1(outlet[i].V[j], outlet[i].r, dir, unit_c.U); - stim::circle c2(outlet[i].V[j + 1], outlet[i].r, dir, unit_c.U); + stim::circle c1(outlet[i].V[j], outlet[i].r * scale, dir, unit_c.U); + stim::circle c2(outlet[i].V[j + 1], outlet[i].r * scale, dir, unit_c.U); cp1 = c1.glpoints(subdivision); cp2 = c2.glpoints(subdivision); @@ -1489,7 +1531,7 @@ namespace stim { } // mark the vertex - void mark_vertex(T scale) { + void mark_vertex(T scale = 1.0f) { glColor3f(0.0f, 0.0f, 0.0f); for (unsigned i = 0; i < num_vertex; i++) { @@ -1500,6 +1542,18 @@ namespace stim { } } + // mark the edge + void mark_edge() { + + glColor3f(0.0f, 1.0f, 0.0f); + for (unsigned i = 0; i < num_edge; i++) { + glRasterPos3f((V[E[i].v[0]] + V[E[i].v[1]])[0]/2, (V[E[i].v[0]] + V[E[i].v[1]])[1] / 2, (V[E[i].v[0]] + V[E[i].v[1]])[2] / 2); + std::stringstream ss; + ss << i; + glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss.str().c_str())); + } + } + // find the nearest vertex of current click position // return true and a value if found inline bool epsilon_vertex(T x, T y, T z, T eps, T scale, unsigned& v) { @@ -2381,7 +2435,7 @@ namespace stim { for (unsigned j = 0; j < num; j++) { if (i != j) { if (inlet[i].V[0][1] == inlet[j].V[0][1]) { - if ((inlet[i].V[1][0] >= inlet[j].V[1][0]) && (fabs(inlet[i].V[1][1]) >= fabs(inlet[j].V[1][1])) && (((inlet[i].V[1][1] - main_feeder[0][1]) * (inlet[j].V[1][1] - main_feeder[0][1])) > 0 ? 1 : 0)) { + if ((inlet[i].V[1][0] >= inlet[j].V[1][0]) && (fabs(inlet[i].V[1][1]) >= fabs(inlet[j].V[1][1])) && (((inlet[i].V[1][1] - main_feeder[0][1]) * (inlet[j].V[1][1] - main_feeder[0][1])) > 0 ? 1 : 0) && inlet[i].V[1][2] == inlet[j].V[1][2]) { inlet_feasibility[i] = false; break; } @@ -2398,7 +2452,7 @@ namespace stim { for (unsigned j = 0; j < num; j++) { if (i != j) { if (outlet[i].V[0][2] == outlet[j].V[0][2]) { - if ((outlet[i].V[1][0] <= outlet[j].V[1][0]) && (fabs(outlet[i].V[1][1]) >= fabs(outlet[j].V[1][1])) && (((outlet[i].V[1][1] - main_feeder[1][1]) * (outlet[j].V[1][1] - main_feeder[1][1])) > 0 ? 1 : 0)) { + if ((outlet[i].V[1][0] <= outlet[j].V[1][0]) && (fabs(outlet[i].V[1][1]) >= fabs(outlet[j].V[1][1])) && (((outlet[i].V[1][1] - main_feeder[1][1]) * (outlet[j].V[1][1] - main_feeder[1][1])) > 0 ? 1 : 0) && outlet[i].V[1][2] == outlet[j].V[1][2]) { outlet_feasibility[i] = false; break; } @@ -2661,6 +2715,11 @@ namespace stim { top = A[i].c[1]; if (A[i].c[1] < bottom) bottom = A[i].c[1]; + // extend the network boundingbox if the additional connections are outside + if (A[i].c[2] > bb.B[2]) + bb.B[2] = A[i].c[2]; + if (A[i].c[2] < bb.A[2]) + bb.A[2] = A[i].c[2]; if (A[i].r > max_radius) max_radius = A[i].r; } diff --git a/main.cu b/main.cu index b476d6f..02a560c 100644 --- a/main.cu +++ b/main.cu @@ -49,7 +49,7 @@ int mods; // special keyboard input std::vector color; // velocity color map std::vector velocity_bar; // velocity bar float length = 40.0f; // cuboid length -float scale = 1.0f; // render scale factor +float scale = 1.0f; // scale factor bool image_stack = false; // flag indicates an image stack been loaded stim::image_stack S; // image stack float binary_threshold = 128; // threshold for binary transformation @@ -57,6 +57,8 @@ float in = 0.0f; // total input volume flow rate float out = 0.0f; float Rt = 0.0f; // total resistance float Qn = 0.0f; // required input total volume flow rate +GLint dlist; // simulation display list +bool undo = false; // delete display list // hard-coded parameters float camera_factor = 1.2f; // start point of the camera as a function of X and Y size @@ -85,7 +87,8 @@ bool simulation = false; // flag indicates simulation mode bool color_bound = false; // flag indicates velocity color map bound bool to_select_pressure = false; // flag indicates having selected a vertex to modify pressure bool mark_index = true; // flag indicates marking the index near the vertex -bool flow_direction = false; // flag indicates rendering glyph for flow velocity field +bool glyph_mode = false; // flag indicates rendering glyph for flow velocity field +bool frame_mode = false; // flag indicates rendering filament framing structrues unsigned pressure_index; // the index of vertex that is clicked unsigned direction_index = -1; // the index of edge that is pointed at unsigned index_index = -1; // the index of the vertex @@ -99,6 +102,7 @@ bool picked_connection = false; // flag indicates picked one connection bool render_new_connection = false; // flag indicates rendering new line connection in trasparency bool redisplay = false; // flag indicates redisplay rendering bool connection_done = false; // flag indicates finishing connections +bool render_flow_rate = false; // flag indicates rendering total volume flow rate unsigned connection_index = -1; // the index of connection that is picked unsigned port_index = 0; // inlet (0) or outlet (1) stim::vec3 tmp_v1, tmp_v2; // temp vertex @@ -112,7 +116,7 @@ bool manufacture = false; // flag indicates manufacture mode // get the network basic information inline void get_background() { - pendant_vertex = flow.get_boundary_vertex(); + pendant_vertex = flow.get_pendant_vertex(); num_edge = flow.edges(); num_vertex = flow.vertices(); @@ -153,9 +157,9 @@ __global__ void binary_transform(unsigned N, T* ptr, F threshold) { if (ix >= N) return; // avoid seg-fault if (ptr[ix] >= threshold) // binary transformation - ptr[ix] = 255; - else ptr[ix] = 0; + else + ptr[ix] = 255; } #endif @@ -252,33 +256,60 @@ void glut_render() { if (!simulation && !build_inlet_outlet || manufacture) { glColor3f(0.0f, 0.0f, 0.0f); - flow.glCylinder0(); + flow.glCylinder0(scale, undo); } - else { - flow.bounding_box(); - // render network - if (!flow_direction) { - flow.glSolidSphere(max_pressure, scale, subdivision); - if (mark_index) - flow.mark_vertex(scale); - //flow.glSolidCone(subdivision); - flow.glSolidCylinder(direction_index, color, scale, subdivision); + else { + flow.bounding_box(); // bounding box + if (num_vertex > 100) { // if the network is big enough (say 100), use display list + if (undo) { // undo rendering list + undo = false; + glDeleteLists(dlist, 1); + } + if (!glIsList(dlist)) { + dlist = glGenLists(1); + glNewList(dlist, GL_COMPILE); + // render network + if (!glyph_mode) { + flow.glSolidSphere(max_pressure, subdivision, scale); + if (mark_index) + flow.mark_vertex(scale); + //flow.glSolidCone(subdivision); + flow.glSolidCylinder(direction_index, color, subdivision, scale); + } + // render glyphs + else + flow.glyph(color, subdivision, scale, frame_mode); + + glEndList(); + } + glCallList(dlist); } - // render glyphs - else - flow.glyph(color, subdivision); - - flow.glSolidCuboid(flow_direction && build_inlet_outlet, subdivision, manufacture, length); - if (render_direction) - flow.glSolidCone(direction_index, scale, subdivision); + else { // small network + // render network + if (!glyph_mode) { + flow.glSolidSphere(max_pressure, subdivision, scale); + if (mark_index) { + flow.mark_vertex(scale); + //flow.mark_edge(); + } + //flow.glSolidCone(subdivision); + flow.glSolidCylinder(direction_index, color, subdivision, scale); + } + // render glyphs + else + flow.glyph(color, subdivision, scale, frame_mode); + } + flow.glSolidCuboid(subdivision, manufacture, length); // render bus source + if (render_direction) // render the flow direction of the vertex pointed + flow.glSolidCone(direction_index, subdivision, scale); } if (build_inlet_outlet) - flow.line_bridge(redisplay, flow_direction); + flow.line_bridge(redisplay); if (manufacture) { - flow.glSolidCuboid(flow_direction, subdivision, manufacture, length); - flow.tube_bridge(subdivision); + flow.glSolidCuboid(subdivision, manufacture, length); + flow.tube_bridge(redisplay, subdivision, scale); } if (picked_connection && render_new_connection) { @@ -422,14 +453,14 @@ void glut_render() { glEnd(); glFlush(); - // pressure bar text + // ratio bar text glColor3f(0.0f, 0.0f, 0.0f); glRasterPos2f(0.0f, vY - border_factor); std::stringstream ss_p; ss_p << "Ratio bar"; glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss_p.str().c_str())); - // pressure range text + // ratio range text float step = vY - 3.0f * border_factor; step /= 10; for (unsigned i = 0; i < 11; i++) { @@ -444,7 +475,8 @@ void glut_render() { } if (build_inlet_outlet) - flow.display_flow_rate(in, out); + if (render_flow_rate) + flow.display_flow_rate(in, out); glutSwapBuffers(); } @@ -456,12 +488,15 @@ void glut_menu(int value) { if (value == 1) { simulation = true; build_inlet_outlet = false; + render_flow_rate = false; manufacture = false; modified_bridge = false; + change_fragment = false; connection_done = false; // first time if (!flow.set) { // only first time simulation called "simulation", ^_^ - flow_initialize(); + get_background(); // get the graph information + flow_initialize(); // initialize flow condition menu_option[0] = "resimulation"; } @@ -469,7 +504,7 @@ void glut_menu(int value) { flow_stable_state(); // main function of solving the linear system flow.print_flow(); - if (!flow_direction) + if (!glyph_mode) glut_set_menu(num, 2); } @@ -496,13 +531,14 @@ void glut_menu(int value) { simulation = false; build_inlet_outlet = false; manufacture = true; - flow_direction = false; // manufacuture mode doesn't need flow direction + glyph_mode = false; // manufacuture mode doesn't need flow direction redisplay = true; } if (value == 4) { simulation = true; build_inlet_outlet = false; + render_flow_rate = false; manufacture = false; adjustment(); // adjust network flow accordingly @@ -541,9 +577,9 @@ void glut_passive_motion(int x, int y) { if (simulation || build_inlet_outlet && !mods) { bool flag = flow.epsilon_edge((float)posX, (float)posY, (float)posZ, eps, direction_index); - if (flag && !flow_direction) + if (flag && !glyph_mode) render_direction = true; - else if (!flag && !flow_direction) { + else if (!flag && !glyph_mode) { if (render_direction) // if the direction is displaying currently, do a short delay Sleep(300); render_direction = false; @@ -654,6 +690,7 @@ void glut_mouse(int button, int state, int x, int y) { fragment_ratio = 0.0f; change_fragment = false; + render_flow_rate = true; flow.modify_synthetic_connection(u, rou, hilbert_curve, height_threshold, in, out, fragment_ratio, default_radius); } // move connections along y-axis @@ -799,7 +836,7 @@ void glut_wheel(int wheel, int direction, int x, int y) { if (!to_select_pressure) { bool flag = flow.epsilon_vertex((float)posX, (float)posY, (float)posZ, eps, scale, pressure_index); - if (flag && simulation) { + if (flag && simulation && !glyph_mode) { float tmp_r; if (direction > 0) { // increase radii tmp_r = flow.get_radius(pressure_index); @@ -812,6 +849,7 @@ void glut_wheel(int wheel, int direction, int x, int y) { tmp_r = default_radius; } flow.set_radius(pressure_index, tmp_r); + undo = true; // undo rendering } else if (!mods) { if (direction > 0) // if it is button 3(up), move closer @@ -839,6 +877,8 @@ void glut_wheel(int wheel, int direction, int x, int y) { else scale = 1.0f; } + undo = true; + redisplay = true; } glutPostRedisplay(); @@ -855,22 +895,44 @@ void glut_keyboard(unsigned char key, int x, int y) { flow.save_network(); break; + // convert network to binary format (.nwt) + case 'c': { + std::vector tmp = stim::parser::split(args.arg(0), '.'); + std::stringstream ss; + ss << tmp[0] << ".nwt"; + std::string filename = ss.str(); + flow.saveNwt(filename); + break; + } + // flow vector field visualization, Glyphs case 'f': - if (flow_direction && !manufacture && (simulation || build_inlet_outlet)) { - flow_direction = false; + if (glyph_mode && !manufacture && (simulation || build_inlet_outlet)) { + glyph_mode = false; + frame_mode = false; redisplay = true; // lines and arrows rendering use the same display list int num = glutGet(GLUT_MENU_NUM_ITEMS); if (num == 1) glut_set_menu(num, 2); } - else if (!flow_direction && !manufacture && (simulation || build_inlet_outlet)) { - flow_direction = true; + else if (!glyph_mode && !manufacture && (simulation || build_inlet_outlet)) { + glyph_mode = true; redisplay = true; int num = glutGet(GLUT_MENU_NUM_ITEMS); if (num == 2) glut_set_menu(num, 1); } + undo = true; + break; + + // filaments around arrows + case 'g': + if (glyph_mode) { + if (frame_mode) + frame_mode = false; + else + frame_mode = true; + } break; // open/close index marks @@ -879,6 +941,7 @@ void glut_keyboard(unsigned char key, int x, int y) { mark_index = false; else mark_index = true; + undo = true; break; // output image stack @@ -947,8 +1010,10 @@ void advertise() { std::cout << "Usage(keyboard): e -> open/close indexing" << std::endl; std::cout << " m -> build synthetic connections(connection mode)/output augmented network as image stack (manufacture mode)" << std::endl; - std::cout << " s -> save network flow profiles in result folder as cvs files" << std::endl; - std::cout << " f -> open/close vector field visualization" << std::endl; + std::cout << " s -> save network flow profiles in profile folder as cvs files" << std::endl; + std::cout << " c -> convert .obj network to .nwt network and stores in main folder" << std::endl; + std::cout << " f -> open/close vector field visualization mode" << std::endl; + std::cout << " g -> render filament frames in vector fiedl visualization mode" << std::endl; std::cout << args.str(); } @@ -966,6 +1031,8 @@ int main(int argc, char* argv[]) { args.add("stack", "load the image stack"); args.add("stackres", "spacing between pixel samples in each dimension(in units/pixel)", "1 1 1", "real value > 0"); args.add("stackdir", "set the directory of the output image stack", "", "any existing directory (ex. /home/name/network)"); + args.add("scale", "scale down rendering fibers"); + args.add("lcc", "extract the largest connected component"); args.parse(argc, argv); // parse the command line @@ -983,6 +1050,8 @@ int main(int argc, char* argv[]) { std::vector tmp = stim::parser::split(args.arg(0), '.'); if ("obj" == tmp[1]) flow.load_obj(args.arg(0)); + else if ("nwt" == tmp[1]) // stim network binary format + flow.loadNwt(args.arg(0)); else if ("swc" == tmp[1]) flow.load_swc(args.arg(0)); else { @@ -990,7 +1059,8 @@ int main(int argc, char* argv[]) { std::exit(1); } } - get_background(); + + // extract the largest connected component // get the units to work on units = args["units"].as_string(); @@ -1016,22 +1086,19 @@ int main(int argc, char* argv[]) { S.load_images(args["stack"].as_string()); // binary transformation #ifdef __CUDACC__ - unsigned N = S.samples(); // number of pixels loaded - unsigned char* d_S; // image stack stored in device + unsigned N = S.samples(); // number of pixels loaded + unsigned char* d_S; // image stack stored in device unsigned char* h_S = (unsigned char*)malloc(N * sizeof(unsigned char)); // image stack stored in host cudaMalloc((void**)&d_S, N * sizeof(unsigned char)); cudaMemcpy(d_S, S.data(), N * sizeof(unsigned char), cudaMemcpyHostToDevice); - unsigned thread = 1024; unsigned block = N / thread + 1; - binary_transform <<>> (N, d_S, binary_threshold); + binary_transform <<>> (N, d_S, binary_threshold); // binaryzation cudaMemcpy(h_S, d_S, N * sizeof(unsigned char), cudaMemcpyDeviceToHost); S.copy(h_S); - - S.save_images("image????.bmp"); #endif } @@ -1044,6 +1111,10 @@ int main(int argc, char* argv[]) { if (args["stackdir"].is_set()) stackdir = args["stackdir"].as_string(); + // get the scale-down factor is provided + if (args["scale"].is_set()) + scale = args["scale"].as_float(); + // glut main loop bb = flow.boundingbox(); glut_initialize(); -- libgit2 0.21.4