From 9e60c6a754f8cefc11dfc2cb509b6274405c4f17 Mon Sep 17 00:00:00 2001 From: Jiaming Guo Date: Fri, 19 May 2017 15:14:08 -0500 Subject: [PATCH] fixed minor bug in edge detection --- flow.h | 323 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------ main.cu | 66 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------- 2 files changed, 328 insertions(+), 61 deletions(-) diff --git a/flow.h b/flow.h index a9bb9ca..ecaac39 100644 --- a/flow.h +++ b/flow.h @@ -331,6 +331,11 @@ namespace stim { return E[tmp_e].r(tmp_v); } + // get the radius of index j of edge i + T get_radius(unsigned i, unsigned j) { + return E[i].r(j); + } + // get the velocity of pendant vertex i T get_velocity(unsigned i) { @@ -495,6 +500,7 @@ namespace stim { order = o; } + // move hilbert curves void move(unsigned i, T *c, direction dir, T dl, int feeder, bool invert) { int cof = (invert) ? -1 : 1; @@ -524,6 +530,7 @@ namespace stim { outlet[i].V.push_back(tmp); } + // form hilbert curves void hilbert_curve(unsigned i, T *c, int order, T dl, int feeder, bool invert, direction dir = DOWN) { if (order == 1) { @@ -727,51 +734,84 @@ namespace stim { // waste? for (unsigned i = 0; i < num_edge; i++) { - for (unsigned j = 0; j < E[i].size(); j++) { - if (j == 0) { // for start vertex - if (P[E[i].v[0]] != 0) { - stim::vec3 new_color; - new_color[0] = (P[E[i].v[0]] / max_pressure) > 0.5f ? 1.0f : 2.0f * P[E[i].v[0]] / max_pressure; // red - new_color[1] = 0.0f; // green - new_color[2] = (P[E[i].v[0]] / max_pressure) > 0.5f ? 1.0f - 2.0f * (P[E[i].v[0]] / max_pressure - 0.5f) : 1.0f; // blue - glColor3f(new_color[0], new_color[1], new_color[2]); - } - } - else if (j == E[i].size() - 1) { // for end vertex - if (P[E[i].v[1]] != 0) { - stim::vec3 new_color; - new_color[0] = (P[E[i].v[1]] / max_pressure) > 0.5f ? 1.0f : 2.0f * P[E[i].v[1]] / max_pressure; // red - new_color[1] = 0.0f; // green - new_color[2] = (P[E[i].v[1]] / max_pressure) > 0.5f ? 1.0f - 2.0f * (P[E[i].v[1]] / max_pressure - 0.5f) : 1.0f; // blue - glColor3f(new_color[0], new_color[1], new_color[2]); - } - } - else - glColor3f(0.5f, 0.5f, 0.5f); // gray point + // draw the starting vertex + if (P[E[i].v[0]] != 0) { + stim::vec3 new_color; + new_color[0] = (P[E[i].v[0]] / max_pressure) > 0.5f ? 1.0f : 2.0f * P[E[i].v[0]] / max_pressure; // red + new_color[1] = 0.0f; // green + new_color[2] = (P[E[i].v[0]] / max_pressure) > 0.5f ? 1.0f - 2.0f * (P[E[i].v[0]] / max_pressure - 0.5f) : 1.0f; // blue + glColor3f(new_color[0], new_color[1], new_color[2]); glPushMatrix(); - glTranslatef(E[i][j][0], E[i][j][1], E[i][j][2]); - glutSolidSphere(get_r(i, j), subdivision, subdivision); + glTranslatef(E[i][0][0], E[i][0][1], E[i][0][2]); + glutSolidSphere(get_r(i, 0), subdivision, subdivision); glPopMatrix(); } + + // draw the ending vertex + if (P[E[i].v[1]] != 0) { + stim::vec3 new_color; + new_color[0] = (P[E[i].v[1]] / max_pressure) > 0.5f ? 1.0f : 2.0f * P[E[i].v[1]] / max_pressure; // red + new_color[1] = 0.0f; // green + new_color[2] = (P[E[i].v[1]] / max_pressure) > 0.5f ? 1.0f - 2.0f * (P[E[i].v[1]] / max_pressure - 0.5f) : 1.0f; // blue + glColor3f(new_color[0], new_color[1], new_color[2]); + + glPushMatrix(); + glTranslatef(E[i][E[i].size() - 1][0], E[i][E[i].size() - 1][1], E[i][E[i].size() - 1][2]); + glutSolidSphere(get_r(i, E[i].size() - 1), subdivision, subdivision); + glPopMatrix(); + } + //for (unsigned j = 0; j < E[i].size(); j++) { + // if (j == 0) { // for start vertex + // if (P[E[i].v[0]] != 0) { + // stim::vec3 new_color; + // new_color[0] = (P[E[i].v[0]] / max_pressure) > 0.5f ? 1.0f : 2.0f * P[E[i].v[0]] / max_pressure; // red + // new_color[1] = 0.0f; // green + // new_color[2] = (P[E[i].v[0]] / max_pressure) > 0.5f ? 1.0f - 2.0f * (P[E[i].v[0]] / max_pressure - 0.5f) : 1.0f; // blue + // glColor3f(new_color[0], new_color[1], new_color[2]); + // } + // } + // else if (j == E[i].size() - 1) { // for end vertex + // if (P[E[i].v[1]] != 0) { + // stim::vec3 new_color; + // new_color[0] = (P[E[i].v[1]] / max_pressure) > 0.5f ? 1.0f : 2.0f * P[E[i].v[1]] / max_pressure; // red + // new_color[1] = 0.0f; // green + // new_color[2] = (P[E[i].v[1]] / max_pressure) > 0.5f ? 1.0f - 2.0f * (P[E[i].v[1]] / max_pressure - 0.5f) : 1.0f; // blue + // glColor3f(new_color[0], new_color[1], new_color[2]); + // } + // } + // else + // glColor3f(0.5f, 0.5f, 0.5f); // gray point + + // glPushMatrix(); + // glTranslatef(E[i][j][0], E[i][j][1], E[i][j][2]); + // glutSolidSphere(get_r(i, j), subdivision, subdivision); + // glPopMatrix(); + //} } } // draw edges as series of cylinders - void glSolidCylinder(GLint subdivision, std::vector color) { + void glSolidCylinder(unsigned index, std::vector color, GLint subdivision) { stim::vec3 tmp_d; stim::vec3 center1; stim::vec3 center2; + stim::circle tmp_c; float r1; float r2; std::vector > cp1(subdivision + 1); std::vector > cp2(subdivision + 1); for (unsigned i = 0; i < num_edge; i++) { // for every edge - glEnable(GL_BLEND); // enable color blend - glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); // set blend function - glDisable(GL_DEPTH_TEST); - glColor4f((float)color[i * 3 + 0] / 255, (float)color[i * 3 + 1] / 255, (float)color[i * 3 + 2] / 255, 0.5f); + if (i == index) { // render in tranparency for direction indication + glEnable(GL_BLEND); // enable color blend + glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); // set blend function + glDisable(GL_DEPTH_TEST); + glColor4f((float)color[i * 3 + 0] / 255, (float)color[i * 3 + 1] / 255, (float)color[i * 3 + 2] / 255, 0.5f); + } + else + glColor3f((float)color[i * 3 + 0] / 255, (float)color[i * 3 + 1] / 255, (float)color[i * 3 + 2] / 255); + for (unsigned j = 0; j < E[i].size() - 1; j++) { // for every point on the edge center1 = E[i][j]; center2 = E[i][j + 1]; @@ -779,8 +819,48 @@ namespace stim { r1 = get_r(i, j); r2 = get_r(i, j + 1); - // calculate the envelope caps - find_envelope(cp1, cp2, center1, center2, r1, r2, subdivision); + //// calculate the envelope caps + //find_envelope(cp1, cp2, center1, center2, r1, r2, subdivision); + 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); + } glBegin(GL_QUAD_STRIP); for (unsigned j = 0; j < cp1.size(); j++) { @@ -789,12 +869,71 @@ namespace stim { } glEnd(); } + if (i == index) { + glDisable(GL_BLEND); + glEnable(GL_DEPTH_TEST); + } + } glFlush(); - glDisable(GL_BLEND); } // draw the flow direction as cone + void glSolidCone(unsigned i, GLint subdivision) { + + stim::vec3 tmp_d; // direction + stim::vec3 center; // cone hat center + stim::vec3 head; // cone hat top + stim::circle tmp_c; + std::vector > cp; + T radius; + + glColor3f(0.0f, 0.0f, 0.0f); + + unsigned index = E[i].size() / 2 - 1; + tmp_d = E[i][index + 1] - E[i][index]; + 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; + if (v[i] > 0) + head = center + tmp_d * 2 * sqrt(3) * radius; + else + head = center - tmp_d * 2 * sqrt(3) * radius; + + 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(); + + // 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]; + // tmp_d = tmp_d.norm(); + // center = (E[i][j + 1] + E[i][j]) / 2; + // tmp_c.rotate(tmp_d); + // radius = (E[i].r(j + 1) + E[i].r(j)) / 2; + // if (v[i] > 0) // if flow flows from j to j+1 + // head = center + tmp_d * 2 * sqrt(3) * radius; + // else + // head = center - tmp_d * 2 * sqrt(3) * radius; + + // 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(); + } void glSolidCone(GLint subdivision) { stim::vec3 tmp_d; // direction @@ -805,27 +944,48 @@ namespace stim { T radius; glColor3f(0.0f, 0.0f, 0.0f); + // draw a cone for every edge to indicate for (unsigned i = 0; i < num_edge; i++) { // for every edge - 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(); - center = (E[i][j + 1] + E[i][j]) / 2; - tmp_c.rotate(tmp_d); - radius = (E[i].r(j + 1) + E[i].r(j)) / 2; - if (v[i] > 0) // if flow flows from j to j+1 - head = center + tmp_d * 2 * sqrt(3) * radius; - else - head = center - tmp_d * 2 * sqrt(3) * radius; - - stim::circle c(center, radius, tmp_d, tmp_c.U); - cp = c.glpoints(subdivision); + unsigned k1 = E[i].size() / 2 - 1; // start and end index + unsigned k2 = E[i].size() / 2; + tmp_d = E[i][k2] - E[i][k1]; + 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; + if (v[i] > 0) // if flow flows from k1 to k2 + head = center + tmp_d * 2 * sqrt(3) * radius; + else + head = center - tmp_d * 2 * sqrt(3) * radius; + 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(); - 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(); + // center = (E[i][j + 1] + E[i][j]) / 2; + // tmp_c.rotate(tmp_d); + // radius = (E[i].r(j + 1) + E[i].r(j)) / 2; + // if (v[i] > 0) // if flow flows from j to j+1 + // head = center + tmp_d * 2 * sqrt(3) * radius; + // else + // head = center - tmp_d * 2 * sqrt(3) * radius; + + // 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(); } @@ -1054,14 +1214,14 @@ namespace stim { // return true and a value if found inline bool epsilon_vertex(T x, T y, T z, T eps, unsigned& v) { - float d = FLT_MAX; // minimum distance between 2 vertices - float tmp_d = 0.0f; // temporary stores distance for loop + T d = FLT_MAX; // minimum distance between 2 vertices + T tmp_d = 0.0f; // temporary stores distance for loop unsigned tmp_i = 0; // temporary stores connection index for loop - stim::vec3 tmp_v; // temporary stores current loop point + stim::vec3 tmp_v; // temporary stores current loop point d = FLT_MAX; // set to max of float number for (unsigned i = 0; i < V.size(); i++) { - tmp_v = stim::vec3(x, y, z); + tmp_v = stim::vec3(x, y, z); tmp_v = tmp_v - V[i]; // calculate a vector between two vertices tmp_d = tmp_v.len(); // calculate length of that vector @@ -1079,6 +1239,61 @@ namespace stim { return false; } + // find the nearest inlet/outlet connection line of current click position + // ab -> line segment, v -> point + // return true and a value if found + inline bool epsilon_edge(T x, T y, T z, T eps, unsigned &idx) { + + T d = FLT_MAX; + T tmp_d; + unsigned tmp_i = 0; + unsigned tmp_j = 0; + stim::vec3 v1; + stim::vec3 v2; + stim::vec3 v3; + stim::vec3 v0 = stim::vec3(x, y, z); + bool online = false; // flag indicates the point is on the line-segment + float a, b; + + for (unsigned i = 0; i < E.size(); i++) { + for (unsigned j = 0; j < E[i].size() - 1; j++) { + v1 = E[i][j + 1] - E[i][j]; // a -> b = ab + v2 = v0 - E[i][j]; // a -> v = av + v3 = v0 - E[i][j + 1]; // b -> v = bv + + tmp_d = v2.dot(v1); // av·ab + + // check the line relative position + a = v2.dot(v1.norm()); + b = v3.dot(v1.norm()); + if (a < v1.len() && b < v1.len()) // if the length of projection fragment is longer than the line-segment + online = true; + else + online = false; + + if (tmp_d <= 0.0 || tmp_d > std::pow(v1.len(), 2) && !online) // projection lies outside the line-segment + continue; + else { + tmp_d = v1.cross(v2).len() / v1.len(); // perpendicular distance of point to segment: |v1 x v2| / |v1| + if (tmp_d < d) { + d = tmp_d; + tmp_i = i; + tmp_j = j; + } + } + } + } + + eps += get_radius(tmp_i, tmp_j); + + if (d < eps) { + idx = tmp_i; + return true; + } + + return false; + } + /// build main feeder connection // set up main feeder and main port of both input and output void set_main_feeder(T border = 400.0f) { diff --git a/main.cu b/main.cu index a122423..caac48a 100644 --- a/main.cu +++ b/main.cu @@ -54,7 +54,7 @@ float zoom_factor = 10.0f; // zooming factor float border_factor = 20.0f; // border float radii_factor = 1.0f; // radii changing factor GLint subdivision = 20; // slices and stacks -float default_radius = 5.0f; // default radii of network vertex +float default_radius = 5.0f; // default radii of network vertex float delta = 0.01f; // small discrepancy float eps = 20.0f; // epsilon threshold float max_pressure = 0.0f; // maximum pressure that the channel can bear @@ -67,10 +67,12 @@ int mouse_y; // window y-coordinate bool LTbutton = false; // true means down while false means up // simulation parameters +bool render_direction = false; // flag indicates rendering flow direction for one edge 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 unsigned pressure_index; // the index of vertex that is clicked +unsigned direction_index = -1; // the index of edge that is pointed at // build inlet/outlet parameters bool build_inlet_outlet = false; // flag indicates building inlets and outlets @@ -92,8 +94,11 @@ inline void get_background() { // set the initial radii flow.init(num_edge, num_vertex); // initialize flow object - for (unsigned i = 0; i < num_edge; i++) - flow.set_r(i, default_radius); + + // if no radius information laoded + if (!flow.get_radius(0, 0)) + for (unsigned i = 0; i < num_edge; i++) + flow.set_r(i, default_radius); } // convert from window coordinates to world coordinates @@ -197,15 +202,18 @@ void glut_render() { glut_modelview(); if (!simulation && !build_inlet_outlet || manufacture) { + glColor3f(0.0f, 0.0f, 0.0f); flow.glCylinder0(); } else { flow.bounding_box(); flow.glSolidSphere(max_pressure, subdivision); flow.mark_vertex(); - flow.glSolidCone(subdivision); - flow.glSolidCylinder(subdivision, color); + //flow.glSolidCone(subdivision); + flow.glSolidCylinder(direction_index, color, subdivision); flow.glSolidCuboid(); + if (render_direction) + flow.glSolidCone(direction_index, subdivision); } if (build_inlet_outlet) { @@ -217,7 +225,6 @@ void glut_render() { flow.tube_bridge(subdivision); } - // render bars // bring up a pressure bar on left if (to_select_pressure) { @@ -416,6 +423,31 @@ void glut_motion(int x, int y) { glutPostRedisplay(); // re-draw the visualization } +// defines passive mouse motion function +void glut_passive_motion(int x, int y) { + + mouse_x = x; + mouse_y = y; + + // check whether the mouse point near to an edge + GLdouble posX, posY, posZ; + window_to_world(posX, posY, posZ); // get the world coordinates + + if (simulation || build_inlet_outlet) { + bool flag = flow.epsilon_edge((float)posX, (float)posY, (float)posZ, eps, direction_index); + if (flag) + render_direction = true; + else { + if (render_direction) // if the direction is displaying currently, do a short delay + Sleep(1000); + render_direction = false; + direction_index = -1; + } + } + + glutPostRedisplay(); // re-draw the visualization +} + // get click window coordinates void glut_mouse(int button, int state, int x, int y) { @@ -552,6 +584,7 @@ void glut_initialize() { glutDisplayFunc(glut_render); glutMouseFunc(glut_mouse); glutMotionFunc(glut_motion); + glutPassiveMotionFunc(glut_passive_motion); glutMouseWheelFunc(glut_wheel); glutKeyboardFunc(glut_keyboard); @@ -565,6 +598,20 @@ void glut_initialize() { cam.LookAt(c[0], c[1], c[2]); } +// output an advertisement for the lab, authors and usage information +void advertise() { + std::cout << std::endl << std::endl; + std::cout << " =======================================================================================" << std::endl; + std::cout << "|Thank you for using the synthetic microvascular model generator for microfluidics tool!|" << std::endl; + std::cout << "|Scalable Tissue Imaging and Modeling (STIM) Lab, University of Houston |" << std::endl; + std::cout << "|Developers: Jiaming Guo, David Mayerich |" << std::endl; + std::cout << "|Source: https://git.stim.ee.uh.edu/Jack/flow3.git |" << std::endl; + std::cout << " =======================================================================================" << std::endl << std::endl; + + std::cout << args.str(); +} + +// main function: parse arguments and initialize GLUT int main(int argc, char* argv[]) { // add arguments @@ -573,12 +620,17 @@ int main(int argc, char* argv[]) { args.add("maxpress", "maximum allowed pressure in g / units / s^2, default 2 is for blood when units = um", "2", "real value > 0"); args.add("viscosity", "set the viscosity of the fluid (in g / units / s), default .00001 is for blood when units = um", ".00001", "real value > 0"); args.add("rou", "set the desity of the fluid (in g / units^3), default 1.06*10^-12 is for blood when units = um", ".00000000000106", "real value > 0"); - args.add("hilbert", "enable hilbert curves connections", "0", "value 1 for enablement"); + args.add("hilbert", "activate hilbert curves connections", "0", "value 1 for enablement"); 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.parse(argc, argv); // parse the command line + if (args["help"].is_set()) { + advertise(); + std::exit(1); + } + // load network if (args["network"].is_set()) { // load network from user std::vector tmp = stim::parser::split(args["network"].as_string(), '.'); -- libgit2 0.21.4