Jiaming Guo / flow3

Commit ce1a6f5e411a92726fc1887f604fc18df15ebe01

Authored by Jiaming Guo
1 parent b710cea4

add functions for incorporating nwt files and prepare for computing largest connected component

Showing 2 changed files with 416 additions and 286 deletions   Show diff stats
Inline Side-by-side
... ... @@ -355,6 +355,24 @@ namespace stim {
355 355 return E[i].r(j);
356 356 }
357 357  
  358 + // get the pendant vertices
  359 + std::vector<unsigned> get_pendant_vertex() {
  360 + std::vector<unsigned> result;
  361 + int count = 0;
  362 +
  363 + for (unsigned i = 0; i < V.size(); i++) { // for every vertex
  364 + for (unsigned j = 0; j < E.size(); j++) { // for every edge
  365 + if (i == E[j].v[0] || i == E[j].v[1]) // check whether current vertex terminates one edge
  366 + count++;
  367 + }
  368 + if (count == 1) // is pendant vertex
  369 + result.push_back(i);
  370 + count = 0; // reset count
  371 + }
  372 +
  373 + return result;
  374 + }
  375 +
358 376 // get the velocity of pendant vertex i
359 377 T get_velocity(unsigned i) {
360 378  
... ... @@ -378,6 +396,12 @@ namespace stim {
378 396 P[i] = value;
379 397 }
380 398  
  399 + // extrct the largest connected component
  400 + void extract_lcc() {
  401 +
  402 +
  403 + }
  404 +
381 405 // solve the linear system to get stable flow state
382 406 void solve_flow(T viscosity) {
383 407  
... ... @@ -385,7 +409,7 @@ namespace stim {
385 409 clear();
386 410  
387 411 // get the pendant vertex indices
388   - pendant_vertex = get_boundary_vertex();
  412 + pendant_vertex = get_pendant_vertex();
389 413  
390 414 // get bounding box
391 415 bb = (*this).boundingbox();
... ... @@ -749,7 +773,7 @@ namespace stim {
749 773 }
750 774  
751 775 // draw solid sphere at every vertex
752   - void glSolidSphere(T max_pressure, T scale, GLint subdivision) {
  776 + void glSolidSphere(T max_pressure, GLint subdivision, T scale = 1.0f) {
753 777  
754 778 // waste?
755 779 for (unsigned i = 0; i < num_edge; i++) {
... ... @@ -808,7 +832,7 @@ namespace stim {
808 832 }
809 833  
810 834 // draw edges as series of cylinders
811   - void glSolidCylinder(unsigned index, std::vector<unsigned char> color, T scale, GLint subdivision) {
  835 + void glSolidCylinder(unsigned index, std::vector<unsigned char> color, GLint subdivision, T scale = 1.0f) {
812 836  
813 837 stim::vec3<float> tmp_d;
814 838 stim::vec3<float> center1;
... ... @@ -895,7 +919,7 @@ namespace stim {
895 919 }
896 920  
897 921 // draw the flow direction as cone, the size of the cone depends on the length of that edge
898   - void glSolidCone(unsigned i, T scale, GLint subdivision) {
  922 + void glSolidCone(unsigned i, GLint subdivision, T scale = 1.0f, T threshold = 0.01f) {
899 923  
900 924 stim::vec3<T> tmp_d; // direction
901 925 stim::vec3<T> center; // cone hat center
... ... @@ -909,27 +933,28 @@ namespace stim {
909 933  
910 934 unsigned index = E[i].size() / 2 - 1;
911 935 tmp_d = E[i][index + 1] - E[i][index];
912   - h = tmp_d.len() / 3.0f; // get the height base by factor 3
  936 + h = tmp_d.len() / 1.5f; // get the height base by factor 3
913 937 tmp_d = tmp_d.norm();
914 938 center = (E[i][index + 1] + E[i][index]) / 2;
915 939 tmp_c.rotate(tmp_d);
916 940 radius = (E[i].r(index + 1) + E[i].r(index)) / 2 * scale;
917 941 radius = (h / sqrt(3) < radius) ? h / sqrt(3) : radius; // update radius
918   - if (v[i] > 0)
  942 + if (v[i] > threshold)
919 943 head = center + tmp_d * h;
920   - else
  944 + else if (v[i] < -threshold)
921 945 head = center - tmp_d * h;
922 946  
923 947 stim::circle<float> c(center, radius, tmp_d, tmp_c.U);
924 948 cp = c.glpoints(subdivision);
925 949  
926   - glBegin(GL_TRIANGLE_FAN);
927   - glVertex3f(head[0], head[1], head[2]);
928   - for (unsigned k = 0; k < cp.size(); k++)
929   - glVertex3f(cp[k][0], cp[k][1], cp[k][2]);
930   - glEnd();
931   - glFlush();
932   -
  950 + if (v[i] > threshold || v[i] < -threshold) {
  951 + glBegin(GL_TRIANGLE_FAN);
  952 + glVertex3f(head[0], head[1], head[2]);
  953 + for (unsigned k = 0; k < cp.size(); k++)
  954 + glVertex3f(cp[k][0], cp[k][1], cp[k][2]);
  955 + glEnd();
  956 + glFlush();
  957 + }
933 958 // draw a cone for every edge to indicate
934 959 //for (unsigned j = 0; j < E[i].size() - 1; j++) { // for every point on current edge
935 960 // tmp_d = E[i][j + 1] - E[i][j];
... ... @@ -953,7 +978,7 @@ namespace stim {
953 978 //}
954 979 //glFlush();
955 980 }
956   - void glSolidCone(GLint subdivision) {
  981 + void glSolidCone(GLint subdivision, T scale = 1.0f, T threhold = 0.01f) {
957 982  
958 983 stim::vec3<T> tmp_d; // direction
959 984 stim::vec3<T> center; // cone hat center
... ... @@ -973,21 +998,22 @@ namespace stim {
973 998 tmp_d = tmp_d.norm();
974 999 center = (E[i][k2] + E[i][k1]) / 2;
975 1000 tmp_c.rotate(tmp_d);
976   - radius = (E[i].r(k2) + E[i].r(k1)) / 2;
  1001 + radius = (E[i].r(k2) + E[i].r(k1)) / 2 * scale;
977 1002 radius = (h / sqrt(3) < radius) ? h / sqrt(3) : radius; // update radius by height base if necessary
978   - if (v[i] > 0) // if flow flows from k1 to k2
  1003 + if (v[i] > threhold) // if flow flows from k1 to k2
979 1004 head = center + tmp_d * h;
980   - else
  1005 + else if(v[i] < -threhold)
981 1006 head = center - tmp_d * h;
982 1007 stim::circle<float> c(center, radius, tmp_d, tmp_c.U);
983 1008 cp = c.glpoints(subdivision);
984 1009  
985   - glBegin(GL_TRIANGLE_FAN);
986   - glVertex3f(head[0], head[1], head[2]);
987   - for (unsigned k = 0; k < cp.size(); k++)
988   - glVertex3f(cp[k][0], cp[k][1], cp[k][2]);
989   - glEnd();
990   -
  1010 + if (v[i] > threhold || v[i] < -threhold) {
  1011 + glBegin(GL_TRIANGLE_FAN);
  1012 + glVertex3f(head[0], head[1], head[2]);
  1013 + for (unsigned k = 0; k < cp.size(); k++)
  1014 + glVertex3f(cp[k][0], cp[k][1], cp[k][2]);
  1015 + glEnd();
  1016 + }
991 1017 //for (unsigned j = 0; j < E[i].size() - 1; j++) { // for every point on current edge
992 1018 // tmp_d = E[i][j + 1] - E[i][j];
993 1019 // tmp_d = tmp_d.norm();
... ... @@ -1013,7 +1039,7 @@ namespace stim {
1013 1039 }
1014 1040  
1015 1041 // draw main feeder as solid cube
1016   - void glSolidCuboid(bool arrow, GLint subdivision, bool manufacture = false, T length = 40.0f, T height = 10.0f) {
  1042 + void glSolidCuboid(GLint subdivision, bool manufacture = false, T length = 40.0f, T height = 10.0f) {
1017 1043  
1018 1044 T width;
1019 1045 stim::vec3<T> L = bb.A; // get the bottom left corner
... ... @@ -1073,51 +1099,10 @@ namespace stim {
1073 1099 glVertex3f(main_feeder[i][0] + length / 2, main_feeder[i][1] - height / 2, main_feeder[i][2] + width / 2);
1074 1100 glEnd();
1075 1101 }
1076   -
1077   - // render total flow direction
1078   - if (arrow) {
1079   - stim::vec3<T> d = main_feeder[1] - main_feeder[0];
1080   - d = d.norm();
1081   - stim::vec3<T> v1, v2, v3;
1082   - stim::vec3<T> center = bb.center();
1083   - stim::vec3<T> size = bb.size();
1084   - stim::circle<float> tmp_c;
1085   - tmp_c.rotate(d);
1086   - std::vector<typename stim::vec3<float> > cp1(subdivision + 1);
1087   - std::vector<typename stim::vec3<float> > cp2(subdivision + 1);
1088   - v2 = center - stim::vec3<T>(0.0f, size[1] + 10.0f, 0.0f);
1089   - v1 = v2 - stim::vec3<T>(30.0f, 0.0f, 0.0f);
1090   - v3 = v2 + stim::vec3<T>(30.0f, 0.0f, 0.0f);
1091   - v2 = v2 + stim::vec3<T>(20.0f, 0.0f, 0.0f);
1092   -
1093   - // draw the total flow direciton indicating arrow
1094   - stim::circle<T> c1(v1, 5.0f / 2, d, tmp_c.U);
1095   - cp1 = c1.glpoints(subdivision);
1096   - stim::circle<T> c2(v2, 5.0f / 2, d, tmp_c.U);
1097   - cp2 = c2.glpoints(subdivision);
1098   -
1099   - glColor3f(0.0f, 1.0f, 0.0f);
1100   - glBegin(GL_QUAD_STRIP);
1101   - for (unsigned k = 0; k < cp1.size(); k++) {
1102   - glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]);
1103   - glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
1104   - }
1105   - glEnd();
1106   -
1107   - // render the cone part
1108   - stim::circle<T> c3(v2, 5.0f, d, tmp_c.U);
1109   - cp2 = c3.glpoints(subdivision);
1110   - glBegin(GL_TRIANGLE_FAN);
1111   - glVertex3f(v3[0], v3[1], v3[2]);
1112   - for (unsigned k = 0; k < cp2.size(); k++)
1113   - glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
1114   - glEnd();
1115   - }
1116   - glFlush();
1117 1102 }
1118 1103  
1119 1104 // draw flow velocity field, glyph
1120   - void glyph(std::vector<unsigned char> color, GLint subdivision, T r = 4.0f) {
  1105 + void glyph(std::vector<unsigned char> color, GLint subdivision, T scale = 1.0f, bool frame = false, T r = 4.0f, T threshold = 0.01f) {
1121 1106  
1122 1107 // v1----v2-->v3
1123 1108 T k = 4.0f; // quartering
... ... @@ -1127,51 +1112,124 @@ namespace stim {
1127 1112 std::vector<typename stim::vec3<float> > cp1(subdivision + 1);
1128 1113 std::vector<typename stim::vec3<float> > cp2(subdivision + 1);
1129 1114  
  1115 + // rendering the arrows
1130 1116 for (unsigned i = 0; i < num_edge; i++) { // for every edge
1131 1117 glColor3f((float)color[i * 3 + 0] / 255.0f, (float)color[i * 3 + 1] / 255.0f, (float)color[i * 3 + 2] / 255.0f);
1132 1118 for (unsigned j = 0; j < E[i].size() - 1; j++) { // for every point on that edge
1133 1119  
1134 1120 // consider the velocity valuence
1135   - if (v[i] > 0) { // positive, from start point to end point
  1121 + if (v[i] > threshold) { // positive, from start point to end point
1136 1122 v1 = E[i][j];
1137 1123 v3 = E[i][j + 1];
1138 1124 }
1139   - else { // negative, from end point to start point
  1125 + else if (v[i] < -threshold) { // negative, from end point to start point
1140 1126 v1 = E[i][j + 1];
1141 1127 v3 = E[i][j];
1142 1128 }
1143   - d = v3 - v1;
1144   - // place the arrow in the middel of one edge
1145   - v2 = v1 + (1.0f / k * 2.0f) * d; // looks like =->=
1146   - v1 = v1 + (1.0f / k) * d;
1147   - v3 = v3 - (1.0f / k) * d;
1148   - d = d.norm();
1149   - tmp_c.rotate(d);
1150   -
1151   - // render the cylinder part
1152   - stim::circle<T> c1(v1, r / 2, d, tmp_c.U);
1153   - cp1 = c1.glpoints(subdivision);
1154   - stim::circle<T> c2(v2, r / 2, d, tmp_c.U);
1155   - cp2 = c2.glpoints(subdivision);
1156   -
1157   - glBegin(GL_QUAD_STRIP);
1158   - for (unsigned k = 0; k < cp1.size(); k++) {
1159   - glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]);
1160   - glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
  1129 +
  1130 + if (v[i] > threshold || v[i] < -threshold) {
  1131 + d = v3 - v1;
  1132 + // place the arrow in the middel of one edge
  1133 + v2 = v1 + (1.0f / k * 2.0f) * d; // looks like =->=
  1134 + v1 = v1 + (1.0f / k) * d;
  1135 + v3 = v3 - (1.0f / k) * d;
  1136 + d = d.norm();
  1137 + tmp_c.rotate(d);
  1138 +
  1139 + // render the cylinder part
  1140 + stim::circle<T> c1(v1, r / 2 * scale, d, tmp_c.U);
  1141 + cp1 = c1.glpoints(subdivision);
  1142 + stim::circle<T> c2(v2, r / 2 * scale, d, tmp_c.U);
  1143 + cp2 = c2.glpoints(subdivision);
  1144 +
  1145 + glBegin(GL_QUAD_STRIP);
  1146 + for (unsigned k = 0; k < cp1.size(); k++) {
  1147 + glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]);
  1148 + glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
  1149 + }
  1150 + glEnd();
  1151 +
  1152 + // render the cone part
  1153 + stim::circle<T> c3(v2, r * scale, d, tmp_c.U);
  1154 + cp2 = c3.glpoints(subdivision);
  1155 + glBegin(GL_TRIANGLE_FAN);
  1156 + glVertex3f(v3[0], v3[1], v3[2]);
  1157 + for (unsigned k = 0; k < cp2.size(); k++)
  1158 + glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
  1159 + glEnd();
1161 1160 }
1162   - glEnd();
  1161 + }
  1162 + }
1163 1163  
1164   - // render the cone part
1165   - stim::circle<T> c3(v2, r, d, tmp_c.U);
1166   - cp2 = c3.glpoints(subdivision);
1167   - glBegin(GL_TRIANGLE_FAN);
1168   - glVertex3f(v3[0], v3[1], v3[2]);
1169   - for (unsigned k = 0; k < cp2.size(); k++)
1170   - glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
1171   - glEnd();
  1164 + // rendering frames
  1165 + if (frame) {
  1166 + frame = false;
  1167 + stim::vec3<float> center1;
  1168 + stim::vec3<float> center2;
  1169 + stim::vec3<float> tmp_d; // flow direction
  1170 + float r1, r2;
  1171 +
  1172 + for (unsigned i = 0; i < num_edge; i++) {
  1173 + for (unsigned j = 0; j < E[i].size() - 1; j++) {
  1174 + center1 = E[i][j];
  1175 + center2 = E[i][j + 1];
  1176 +
  1177 + r1 = get_r(i, j) * scale;
  1178 + r2 = get_r(i, j + 1) * scale;
  1179 +
  1180 + subdivision = 5; // rough frames
  1181 +
  1182 + if (j == 0) {
  1183 + if (E[i].size() == 2)
  1184 + find_envelope(cp1, cp2, center1, center2, r1, r2, subdivision);
  1185 + else {
  1186 + tmp_d = center2 - center1;
  1187 + tmp_d = tmp_d.norm();
  1188 + tmp_c.rotate(tmp_d);
  1189 + stim::circle<float> c1(center1, r1, tmp_d, tmp_c.U);
  1190 + cp1 = c1.glpoints(subdivision);
  1191 + tmp_d = (E[i][j + 2] - center2) + (center2 - center1);
  1192 + tmp_d = tmp_d.norm();
  1193 + tmp_c.rotate(tmp_d);
  1194 + stim::circle<float> c2(center2, r2, tmp_d, tmp_c.U);
  1195 + cp2 = c2.glpoints(subdivision);
  1196 + }
  1197 + }
  1198 + else if (j == E[i].size() - 2) {
  1199 + tmp_d = (center2 - center1) + (center1 - E[i][j - 1]);
  1200 + tmp_d = tmp_d.norm();
  1201 + tmp_c.rotate(tmp_d);
  1202 + stim::circle<float> c1(center1, r1, tmp_d, tmp_c.U);
  1203 + cp1 = c1.glpoints(subdivision);
  1204 + tmp_d = center2 - center1;
  1205 + tmp_d = tmp_d.norm();
  1206 + tmp_c.rotate(tmp_d);
  1207 + stim::circle<float> c2(center2, r2, tmp_d, tmp_c.U);
  1208 + cp2 = c2.glpoints(subdivision);
  1209 + }
  1210 + else {
  1211 + tmp_d = (center2 - center1) + (center1 - E[i][j - 1]);
  1212 + tmp_d = tmp_d.norm();
  1213 + tmp_c.rotate(tmp_d);
  1214 + stim::circle<float> c1(center1, r1, tmp_d, tmp_c.U);
  1215 + cp1 = c1.glpoints(subdivision);
  1216 + tmp_d = (E[i][j + 2] - center2) + (center2 - center1);
  1217 + tmp_d = tmp_d.norm();
  1218 + tmp_c.rotate(tmp_d);
  1219 + stim::circle<float> c2(center2, r2, tmp_d, tmp_c.U);
  1220 + cp2 = c2.glpoints(subdivision);
  1221 + }
  1222 +
  1223 + glColor3f(140/255.0f, 81/255.0f, 10/255.0f);
  1224 + glBegin(GL_LINES);
  1225 + for (unsigned k = 0; k < cp1.size(); k++) {
  1226 + glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]);
  1227 + glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
  1228 + }
  1229 + glEnd();
  1230 + }
1172 1231 }
1173 1232 }
1174   - glFlush();
1175 1233 }
1176 1234  
1177 1235 // display the total volume flow rate
... ... @@ -1204,145 +1262,121 @@ namespace stim {
1204 1262 }
1205 1263  
1206 1264 // draw the bridge as lines or arrows
1207   - void line_bridge(bool &redisplay, bool arrow, T r = 4.0f) {
  1265 + void line_bridge(bool &redisplay, T r = 4.0f) {
1208 1266  
1209   - if (redisplay)
  1267 + if (redisplay) { // check to see whether the display list needs to be updated
1210 1268 glDeleteLists(dlist, 1);
1211   - redisplay = false;
  1269 + redisplay = false;
  1270 + }
1212 1271  
1213 1272 if (!glIsList(dlist)) {
1214 1273 dlist = glGenLists(1);
1215 1274 glNewList(dlist, GL_COMPILE);
1216 1275  
1217   - // render flow direction arrows
1218   - if (arrow) {
1219   - // v1----v2-->v3
1220   - T k = 4.0f; // quartering
1221   - stim::vec3<T> v1, v2, v3; // three point
1222   - stim::vec3<T> d; // direction vector
1223   - stim::circle<float> tmp_c;
1224   - std::vector<typename stim::vec3<float> > cp1(subdivision + 1);
1225   - std::vector<typename stim::vec3<float> > cp2(subdivision + 1);
1226   -
1227   - // inlet, right-going
1228   - for (unsigned i = 0; i < inlet.size(); i++) {
1229   - if (inlet_feasibility[i])
1230   - glColor3f(0.0f, 0.0f, 0.0f); // feasible -> black
1231   - else
1232   - glColor3f(1.0f, 0.0f, 0.0f); // nonfeasible -> red
1233   - for (unsigned j = 0; j < inlet[i].V.size() - 1; j++) {
1234   - v1 = inlet[i].V[j];
1235   - v3 = inlet[i].V[j + 1];
1236   - d = v3 - v1;
1237   - // place the arrow in the middel of one edge
1238   - v2 = v1 + (1.0f / k * 2.0f) * d; // looks like =->=
1239   - v1 = v1 + (1.0f / k) * d;
1240   - v3 = v3 - (1.0f / k) * d;
1241   - d = d.norm();
1242   - tmp_c.rotate(d);
1243   -
1244   - // render the cylinder part
1245   - stim::circle<T> c1(v1, r / 2, d, tmp_c.U);
1246   - cp1 = c1.glpoints(subdivision);
1247   - stim::circle<T> c2(v2, r / 2, d, tmp_c.U);
1248   - cp2 = c2.glpoints(subdivision);
1249   -
1250   - glBegin(GL_QUAD_STRIP);
1251   - for (unsigned k = 0; k < cp1.size(); k++) {
1252   - glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]);
1253   - glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
1254   - }
1255   - glEnd();
1256   -
1257   - // render the cone part
1258   - stim::circle<T> c3(v2, r, d, tmp_c.U);
1259   - cp2 = c3.glpoints(subdivision);
1260   - glBegin(GL_TRIANGLE_FAN);
1261   - glVertex3f(v3[0], v3[1], v3[2]);
1262   - for (unsigned k = 0; k < cp2.size(); k++)
1263   - glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
1264   - glEnd();
1265   - }
1266   - }
1267   -
1268   - // outlet, right-going
1269   - for (unsigned i = 0; i < outlet.size(); i++) {
1270   - if (outlet_feasibility[i])
1271   - glColor3f(0.0f, 0.0f, 0.0f); // feasible -> black
1272   - else
1273   - glColor3f(1.0f, 0.0f, 0.0f); // nonfeasible -> red
1274   - for (unsigned j = 0; j < outlet[i].V.size() - 1; j++) {
1275   - v1 = outlet[i].V[j + 1];
1276   - v3 = outlet[i].V[j];
1277   - d = v3 - v1;
1278   - // place the arrow in the middel of one edge
1279   - v2 = v1 + (1.0f / k * 2.0f) * d; // looks like =->=
1280   - v1 = v1 + (1.0f / k) * d;
1281   - v3 = v3 - (1.0f / k) * d;
1282   - d = d.norm();
1283   - tmp_c.rotate(d);
1284   -
1285   - // render the cylinder part
1286   - stim::circle<T> c1(v1, r / 2, d, tmp_c.U);
1287   - cp1 = c1.glpoints(subdivision);
1288   - stim::circle<T> c2(v2, r / 2, d, tmp_c.U);
1289   - cp2 = c2.glpoints(subdivision);
1290   -
1291   - glBegin(GL_QUAD_STRIP);
1292   - for (unsigned k = 0; k < cp1.size(); k++) {
1293   - glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]);
1294   - glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
1295   - }
1296   - glEnd();
1297   -
1298   - // render the cone part
1299   - stim::circle<T> c3(v2, r, d, tmp_c.U);
1300   - cp2 = c3.glpoints(subdivision);
1301   - glBegin(GL_TRIANGLE_FAN);
1302   - glVertex3f(v3[0], v3[1], v3[2]);
1303   - for (unsigned k = 0; k < cp2.size(); k++)
1304   - glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
1305   - glEnd();
1306   - }
1307   - }
1308   -
1309   - // render transparent lines as indexing
1310   - glEnable(GL_BLEND); // enable color blend
1311   - glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); // set blend function
1312   - glDisable(GL_DEPTH_TEST);
1313   - glLineWidth(5);
1314   - for (unsigned i = 0; i < inlet.size(); i++) {
1315   - if (inlet_feasibility[i])
1316   - glColor4f(0.0f, 0.0f, 0.0f, 0.2f); // feasible -> black
1317   - else
1318   - glColor4f(1.0f, 0.0f, 0.0f, 0.2f); // nonfeasible -> red
1319   -
1320   - glBegin(GL_LINE_STRIP);
1321   - for (unsigned j = 0; j < inlet[i].V.size(); j++)
1322   - glVertex3f(inlet[i].V[j][0], inlet[i].V[j][1], inlet[i].V[j][2]);
1323   - glEnd();
1324   - }
1325   - for (unsigned i = 0; i < outlet.size(); i++) {
1326   - if (outlet_feasibility[i])
1327   - glColor4f(0.0f, 0.0f, 0.0f, 0.2f); // feasible -> black
1328   - else
1329   - glColor4f(1.0f, 0.0f, 0.0f, 0.2f); // nonfeasible -> red
1330   - glBegin(GL_LINE_STRIP);
1331   - for (unsigned j = 0; j < outlet[i].V.size(); j++)
1332   - glVertex3f(outlet[i].V[j][0], outlet[i].V[j][1], outlet[i].V[j][2]);
1333   - glEnd();
1334   - }
1335   - glDisable(GL_BLEND);
1336   - glEnable(GL_DEPTH_TEST);
1337   - }
1338   - // render connection lines
1339   - else {
  1276 + //// render flow direction arrows
  1277 + //if (arrow) {
  1278 + // // v1----v2-->v3
  1279 + // T k = 4.0f; // quartering
  1280 + // stim::vec3<T> v1, v2, v3; // three point
  1281 + // stim::vec3<T> d; // direction vector
  1282 + // stim::circle<float> tmp_c;
  1283 + // std::vector<typename stim::vec3<float> > cp1(subdivision + 1);
  1284 + // std::vector<typename stim::vec3<float> > cp2(subdivision + 1);
  1285 +
  1286 + // // inlet, right-going
  1287 + // for (unsigned i = 0; i < inlet.size(); i++) {
  1288 + // if (inlet_feasibility[i])
  1289 + // glColor3f(0.0f, 0.0f, 0.0f); // feasible -> black
  1290 + // else
  1291 + // glColor3f(1.0f, 0.0f, 0.0f); // nonfeasible -> red
  1292 + // for (unsigned j = 0; j < inlet[i].V.size() - 1; j++) {
  1293 + // v1 = inlet[i].V[j];
  1294 + // v3 = inlet[i].V[j + 1];
  1295 + // d = v3 - v1;
  1296 + // // place the arrow in the middel of one edge
  1297 + // v2 = v1 + (1.0f / k * 2.0f) * d; // looks like =->=
  1298 + // v1 = v1 + (1.0f / k) * d;
  1299 + // v3 = v3 - (1.0f / k) * d;
  1300 + // d = d.norm();
  1301 + // tmp_c.rotate(d);
  1302 +
  1303 + // // render the cylinder part
  1304 + // stim::circle<T> c1(v1, r / 2, d, tmp_c.U);
  1305 + // cp1 = c1.glpoints(subdivision);
  1306 + // stim::circle<T> c2(v2, r / 2, d, tmp_c.U);
  1307 + // cp2 = c2.glpoints(subdivision);
  1308 +
  1309 + // glBegin(GL_QUAD_STRIP);
  1310 + // for (unsigned k = 0; k < cp1.size(); k++) {
  1311 + // glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]);
  1312 + // glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
  1313 + // }
  1314 + // glEnd();
  1315 +
  1316 + // // render the cone part
  1317 + // stim::circle<T> c3(v2, r, d, tmp_c.U);
  1318 + // cp2 = c3.glpoints(subdivision);
  1319 + // glBegin(GL_TRIANGLE_FAN);
  1320 + // glVertex3f(v3[0], v3[1], v3[2]);
  1321 + // for (unsigned k = 0; k < cp2.size(); k++)
  1322 + // glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
  1323 + // glEnd();
  1324 + // }
  1325 + // }
  1326 +
  1327 + // // outlet, right-going
  1328 + // for (unsigned i = 0; i < outlet.size(); i++) {
  1329 + // if (outlet_feasibility[i])
  1330 + // glColor3f(0.0f, 0.0f, 0.0f); // feasible -> black
  1331 + // else
  1332 + // glColor3f(1.0f, 0.0f, 0.0f); // nonfeasible -> red
  1333 + // for (unsigned j = 0; j < outlet[i].V.size() - 1; j++) {
  1334 + // v1 = outlet[i].V[j + 1];
  1335 + // v3 = outlet[i].V[j];
  1336 + // d = v3 - v1;
  1337 + // // place the arrow in the middel of one edge
  1338 + // v2 = v1 + (1.0f / k * 2.0f) * d; // looks like =->=
  1339 + // v1 = v1 + (1.0f / k) * d;
  1340 + // v3 = v3 - (1.0f / k) * d;
  1341 + // d = d.norm();
  1342 + // tmp_c.rotate(d);
  1343 +
  1344 + // // render the cylinder part
  1345 + // stim::circle<T> c1(v1, r / 2, d, tmp_c.U);
  1346 + // cp1 = c1.glpoints(subdivision);
  1347 + // stim::circle<T> c2(v2, r / 2, d, tmp_c.U);
  1348 + // cp2 = c2.glpoints(subdivision);
  1349 +
  1350 + // glBegin(GL_QUAD_STRIP);
  1351 + // for (unsigned k = 0; k < cp1.size(); k++) {
  1352 + // glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]);
  1353 + // glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
  1354 + // }
  1355 + // glEnd();
  1356 +
  1357 + // // render the cone part
  1358 + // stim::circle<T> c3(v2, r, d, tmp_c.U);
  1359 + // cp2 = c3.glpoints(subdivision);
  1360 + // glBegin(GL_TRIANGLE_FAN);
  1361 + // glVertex3f(v3[0], v3[1], v3[2]);
  1362 + // for (unsigned k = 0; k < cp2.size(); k++)
  1363 + // glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
  1364 + // glEnd();
  1365 + // }
  1366 + // }
  1367 +
  1368 + // // render transparent lines as indexing
  1369 + // glEnable(GL_BLEND); // enable color blend
  1370 + // glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); // set blend function
  1371 + // glDisable(GL_DEPTH_TEST);
1340 1372 glLineWidth(5);
1341 1373 for (unsigned i = 0; i < inlet.size(); i++) {
1342 1374 if (inlet_feasibility[i])
1343   - glColor3f(0.0f, 0.0f, 0.0f); // feasible -> black
  1375 + glColor3f(0.0f, 0.0f, 0.0f);
  1376 + // glColor4f(0.0f, 0.0f, 0.0f, 0.2f); // feasible -> black
1344 1377 else
1345   - glColor3f(1.0f, 0.0f, 0.0f); // nonfeasible -> red
  1378 + glColor3f(1.0f, 0.0f, 0.0f);
  1379 + // glColor4f(1.0f, 0.0f, 0.0f, 0.2f); // nonfeasible -> red
1346 1380  
1347 1381 glBegin(GL_LINE_STRIP);
1348 1382 for (unsigned j = 0; j < inlet[i].V.size(); j++)
... ... @@ -1351,23 +1385,31 @@ namespace stim {
1351 1385 }
1352 1386 for (unsigned i = 0; i < outlet.size(); i++) {
1353 1387 if (outlet_feasibility[i])
1354   - glColor3f(0.0f, 0.0f, 0.0f); // feasible -> black
  1388 + glColor3f(0.0f, 0.0f, 0.0f);
  1389 + // glColor4f(0.0f, 0.0f, 0.0f, 0.2f); // feasible -> black
1355 1390 else
1356   - glColor3f(1.0f, 0.0f, 0.0f); // nonfeasible -> red
  1391 + glColor3f(1.0f, 0.0f, 0.0f);
  1392 + // glColor4f(1.0f, 0.0f, 0.0f, 0.2f); // nonfeasible -> red
1357 1393 glBegin(GL_LINE_STRIP);
1358 1394 for (unsigned j = 0; j < outlet[i].V.size(); j++)
1359 1395 glVertex3f(outlet[i].V[j][0], outlet[i].V[j][1], outlet[i].V[j][2]);
1360 1396 glEnd();
1361 1397 }
1362   - }
1363   - glFlush();
  1398 + // glDisable(GL_BLEND);
  1399 + // glEnable(GL_DEPTH_TEST);
  1400 + //}
1364 1401 glEndList();
1365 1402 }
1366 1403 glCallList(dlist);
1367 1404 }
1368 1405  
1369 1406 // draw the bridge as tubes
1370   - void tube_bridge(T subdivision, T radius = 5.0f) {
  1407 + void tube_bridge(bool &redisplay, T subdivision, T scale = 1.0f, T radius = 5.0f) {
  1408 +
  1409 + if (redisplay) {
  1410 + glDeleteLists(dlist + 1, 1);
  1411 + redisplay = false;
  1412 + }
1371 1413  
1372 1414 if (!glIsList(dlist + 1)) {
1373 1415 glNewList(dlist + 1, GL_COMPILE);
... ... @@ -1383,7 +1425,7 @@ namespace stim {
1383 1425 for (unsigned j = 1; j < inlet[i].V.size() - 1; j++) {
1384 1426 glPushMatrix();
1385 1427 glTranslatef(inlet[i].V[j][0], inlet[i].V[j][1], inlet[i].V[j][2]);
1386   - glutSolidSphere(radius, subdivision, subdivision);
  1428 + glutSolidSphere(radius * scale, subdivision, subdivision);
1387 1429 glPopMatrix();
1388 1430 }
1389 1431 // render edge as cylinder
... ... @@ -1391,8 +1433,8 @@ namespace stim {
1391 1433 dir = inlet[i].V[j] - inlet[i].V[j + 1];
1392 1434 dir = dir.norm();
1393 1435 unit_c.rotate(dir);
1394   - stim::circle<T> c1(inlet[i].V[j], inlet[i].r, dir, unit_c.U);
1395   - stim::circle<T> c2(inlet[i].V[j + 1], inlet[i].r, dir, unit_c.U);
  1436 + stim::circle<T> c1(inlet[i].V[j], inlet[i].r * scale, dir, unit_c.U);
  1437 + stim::circle<T> c2(inlet[i].V[j + 1], inlet[i].r * scale, dir, unit_c.U);
1396 1438 cp1 = c1.glpoints(subdivision);
1397 1439 cp2 = c2.glpoints(subdivision);
1398 1440  
... ... @@ -1410,7 +1452,7 @@ namespace stim {
1410 1452 for (unsigned j = 1; j < outlet[i].V.size() - 1; j++) {
1411 1453 glPushMatrix();
1412 1454 glTranslatef(outlet[i].V[j][0], outlet[i].V[j][1], outlet[i].V[j][2]);
1413   - glutSolidSphere(radius, subdivision, subdivision);
  1455 + glutSolidSphere(radius * scale, subdivision, subdivision);
1414 1456 glPopMatrix();
1415 1457 }
1416 1458 // render edge as cylinder
... ... @@ -1418,8 +1460,8 @@ namespace stim {
1418 1460 dir = outlet[i].V[j] - outlet[i].V[j + 1];
1419 1461 dir = dir.norm();
1420 1462 unit_c.rotate(dir);
1421   - stim::circle<T> c1(outlet[i].V[j], outlet[i].r, dir, unit_c.U);
1422   - stim::circle<T> c2(outlet[i].V[j + 1], outlet[i].r, dir, unit_c.U);
  1463 + stim::circle<T> c1(outlet[i].V[j], outlet[i].r * scale, dir, unit_c.U);
  1464 + stim::circle<T> c2(outlet[i].V[j + 1], outlet[i].r * scale, dir, unit_c.U);
1423 1465 cp1 = c1.glpoints(subdivision);
1424 1466 cp2 = c2.glpoints(subdivision);
1425 1467  
... ... @@ -1489,7 +1531,7 @@ namespace stim {
1489 1531 }
1490 1532  
1491 1533 // mark the vertex
1492   - void mark_vertex(T scale) {
  1534 + void mark_vertex(T scale = 1.0f) {
1493 1535  
1494 1536 glColor3f(0.0f, 0.0f, 0.0f);
1495 1537 for (unsigned i = 0; i < num_vertex; i++) {
... ... @@ -1500,6 +1542,18 @@ namespace stim {
1500 1542 }
1501 1543 }
1502 1544  
  1545 + // mark the edge
  1546 + void mark_edge() {
  1547 +
  1548 + glColor3f(0.0f, 1.0f, 0.0f);
  1549 + for (unsigned i = 0; i < num_edge; i++) {
  1550 + 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);
  1551 + std::stringstream ss;
  1552 + ss << i;
  1553 + glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss.str().c_str()));
  1554 + }
  1555 + }
  1556 +
1503 1557 // find the nearest vertex of current click position
1504 1558 // return true and a value if found
1505 1559 inline bool epsilon_vertex(T x, T y, T z, T eps, T scale, unsigned& v) {
... ... @@ -2381,7 +2435,7 @@ namespace stim {
2381 2435 for (unsigned j = 0; j < num; j++) {
2382 2436 if (i != j) {
2383 2437 if (inlet[i].V[0][1] == inlet[j].V[0][1]) {
2384   - 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)) {
  2438 + 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]) {
2385 2439 inlet_feasibility[i] = false;
2386 2440 break;
2387 2441 }
... ... @@ -2398,7 +2452,7 @@ namespace stim {
2398 2452 for (unsigned j = 0; j < num; j++) {
2399 2453 if (i != j) {
2400 2454 if (outlet[i].V[0][2] == outlet[j].V[0][2]) {
2401   - 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)) {
  2455 + 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]) {
2402 2456 outlet_feasibility[i] = false;
2403 2457 break;
2404 2458 }
... ... @@ -2661,6 +2715,11 @@ namespace stim {
2661 2715 top = A[i].c[1];
2662 2716 if (A[i].c[1] < bottom)
2663 2717 bottom = A[i].c[1];
  2718 + // extend the network boundingbox if the additional connections are outside
  2719 + if (A[i].c[2] > bb.B[2])
  2720 + bb.B[2] = A[i].c[2];
  2721 + if (A[i].c[2] < bb.A[2])
  2722 + bb.A[2] = A[i].c[2];
2664 2723 if (A[i].r > max_radius)
2665 2724 max_radius = A[i].r;
2666 2725 }
... ...
... ... @@ -49,7 +49,7 @@ int mods; // special keyboard input
49 49 std::vector<unsigned char> color; // velocity color map
50 50 std::vector<int> velocity_bar; // velocity bar
51 51 float length = 40.0f; // cuboid length
52   -float scale = 1.0f; // render scale factor
  52 +float scale = 1.0f; // scale factor
53 53 bool image_stack = false; // flag indicates an image stack been loaded
54 54 stim::image_stack<unsigned char, float> S; // image stack
55 55 float binary_threshold = 128; // threshold for binary transformation
... ... @@ -57,6 +57,8 @@ float in = 0.0f; // total input volume flow rate
57 57 float out = 0.0f;
58 58 float Rt = 0.0f; // total resistance
59 59 float Qn = 0.0f; // required input total volume flow rate
  60 +GLint dlist; // simulation display list
  61 +bool undo = false; // delete display list
60 62  
61 63 // hard-coded parameters
62 64 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
85 87 bool color_bound = false; // flag indicates velocity color map bound
86 88 bool to_select_pressure = false; // flag indicates having selected a vertex to modify pressure
87 89 bool mark_index = true; // flag indicates marking the index near the vertex
88   -bool flow_direction = false; // flag indicates rendering glyph for flow velocity field
  90 +bool glyph_mode = false; // flag indicates rendering glyph for flow velocity field
  91 +bool frame_mode = false; // flag indicates rendering filament framing structrues
89 92 unsigned pressure_index; // the index of vertex that is clicked
90 93 unsigned direction_index = -1; // the index of edge that is pointed at
91 94 unsigned index_index = -1; // the index of the vertex
... ... @@ -99,6 +102,7 @@ bool picked_connection = false; // flag indicates picked one connection
99 102 bool render_new_connection = false; // flag indicates rendering new line connection in trasparency
100 103 bool redisplay = false; // flag indicates redisplay rendering
101 104 bool connection_done = false; // flag indicates finishing connections
  105 +bool render_flow_rate = false; // flag indicates rendering total volume flow rate
102 106 unsigned connection_index = -1; // the index of connection that is picked
103 107 unsigned port_index = 0; // inlet (0) or outlet (1)
104 108 stim::vec3<float> tmp_v1, tmp_v2; // temp vertex
... ... @@ -112,7 +116,7 @@ bool manufacture = false; // flag indicates manufacture mode
112 116 // get the network basic information
113 117 inline void get_background() {
114 118  
115   - pendant_vertex = flow.get_boundary_vertex();
  119 + pendant_vertex = flow.get_pendant_vertex();
116 120 num_edge = flow.edges();
117 121 num_vertex = flow.vertices();
118 122  
... ... @@ -153,9 +157,9 @@ __global__ void binary_transform(unsigned N, T* ptr, F threshold) {
153 157 if (ix >= N) return; // avoid seg-fault
154 158  
155 159 if (ptr[ix] >= threshold) // binary transformation
156   - ptr[ix] = 255;
157   - else
158 160 ptr[ix] = 0;
  161 + else
  162 + ptr[ix] = 255;
159 163 }
160 164 #endif
161 165  
... ... @@ -252,33 +256,60 @@ void glut_render() {
252 256  
253 257 if (!simulation && !build_inlet_outlet || manufacture) {
254 258 glColor3f(0.0f, 0.0f, 0.0f);
255   - flow.glCylinder0();
  259 + flow.glCylinder0(scale, undo);
256 260 }
257   - else {
258   - flow.bounding_box();
259   - // render network
260   - if (!flow_direction) {
261   - flow.glSolidSphere(max_pressure, scale, subdivision);
262   - if (mark_index)
263   - flow.mark_vertex(scale);
264   - //flow.glSolidCone(subdivision);
265   - flow.glSolidCylinder(direction_index, color, scale, subdivision);
  261 + else {
  262 + flow.bounding_box(); // bounding box
  263 + if (num_vertex > 100) { // if the network is big enough (say 100), use display list
  264 + if (undo) { // undo rendering list
  265 + undo = false;
  266 + glDeleteLists(dlist, 1);
  267 + }
  268 + if (!glIsList(dlist)) {
  269 + dlist = glGenLists(1);
  270 + glNewList(dlist, GL_COMPILE);
  271 + // render network
  272 + if (!glyph_mode) {
  273 + flow.glSolidSphere(max_pressure, subdivision, scale);
  274 + if (mark_index)
  275 + flow.mark_vertex(scale);
  276 + //flow.glSolidCone(subdivision);
  277 + flow.glSolidCylinder(direction_index, color, subdivision, scale);
  278 + }
  279 + // render glyphs
  280 + else
  281 + flow.glyph(color, subdivision, scale, frame_mode);
  282 +
  283 + glEndList();
  284 + }
  285 + glCallList(dlist);
266 286 }
267   - // render glyphs
268   - else
269   - flow.glyph(color, subdivision);
270   -
271   - flow.glSolidCuboid(flow_direction && build_inlet_outlet, subdivision, manufacture, length);
272   - if (render_direction)
273   - flow.glSolidCone(direction_index, scale, subdivision);
  287 + else { // small network
  288 + // render network
  289 + if (!glyph_mode) {
  290 + flow.glSolidSphere(max_pressure, subdivision, scale);
  291 + if (mark_index) {
  292 + flow.mark_vertex(scale);
  293 + //flow.mark_edge();
  294 + }
  295 + //flow.glSolidCone(subdivision);
  296 + flow.glSolidCylinder(direction_index, color, subdivision, scale);
  297 + }
  298 + // render glyphs
  299 + else
  300 + flow.glyph(color, subdivision, scale, frame_mode);
  301 + }
  302 + flow.glSolidCuboid(subdivision, manufacture, length); // render bus source
  303 + if (render_direction) // render the flow direction of the vertex pointed
  304 + flow.glSolidCone(direction_index, subdivision, scale);
274 305 }
275 306  
276 307 if (build_inlet_outlet)
277   - flow.line_bridge(redisplay, flow_direction);
  308 + flow.line_bridge(redisplay);
278 309  
279 310 if (manufacture) {
280   - flow.glSolidCuboid(flow_direction, subdivision, manufacture, length);
281   - flow.tube_bridge(subdivision);
  311 + flow.glSolidCuboid(subdivision, manufacture, length);
  312 + flow.tube_bridge(redisplay, subdivision, scale);
282 313 }
283 314  
284 315 if (picked_connection && render_new_connection) {
... ... @@ -422,14 +453,14 @@ void glut_render() {
422 453 glEnd();
423 454 glFlush();
424 455  
425   - // pressure bar text
  456 + // ratio bar text
426 457 glColor3f(0.0f, 0.0f, 0.0f);
427 458 glRasterPos2f(0.0f, vY - border_factor);
428 459 std::stringstream ss_p;
429 460 ss_p << "Ratio bar";
430 461 glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss_p.str().c_str()));
431 462  
432   - // pressure range text
  463 + // ratio range text
433 464 float step = vY - 3.0f * border_factor;
434 465 step /= 10;
435 466 for (unsigned i = 0; i < 11; i++) {
... ... @@ -444,7 +475,8 @@ void glut_render() {
444 475 }
445 476  
446 477 if (build_inlet_outlet)
447   - flow.display_flow_rate(in, out);
  478 + if (render_flow_rate)
  479 + flow.display_flow_rate(in, out);
448 480  
449 481 glutSwapBuffers();
450 482 }
... ... @@ -456,12 +488,15 @@ void glut_menu(int value) {
456 488 if (value == 1) {
457 489 simulation = true;
458 490 build_inlet_outlet = false;
  491 + render_flow_rate = false;
459 492 manufacture = false;
460 493 modified_bridge = false;
  494 + change_fragment = false;
461 495 connection_done = false;
462 496 // first time
463 497 if (!flow.set) { // only first time simulation called "simulation", ^_^
464   - flow_initialize();
  498 + get_background(); // get the graph information
  499 + flow_initialize(); // initialize flow condition
465 500 menu_option[0] = "resimulation";
466 501 }
467 502  
... ... @@ -469,7 +504,7 @@ void glut_menu(int value) {
469 504 flow_stable_state(); // main function of solving the linear system
470 505 flow.print_flow();
471 506  
472   - if (!flow_direction)
  507 + if (!glyph_mode)
473 508 glut_set_menu(num, 2);
474 509 }
475 510  
... ... @@ -496,13 +531,14 @@ void glut_menu(int value) {
496 531 simulation = false;
497 532 build_inlet_outlet = false;
498 533 manufacture = true;
499   - flow_direction = false; // manufacuture mode doesn't need flow direction
  534 + glyph_mode = false; // manufacuture mode doesn't need flow direction
500 535 redisplay = true;
501 536 }
502 537  
503 538 if (value == 4) {
504 539 simulation = true;
505 540 build_inlet_outlet = false;
  541 + render_flow_rate = false;
506 542 manufacture = false;
507 543  
508 544 adjustment(); // adjust network flow accordingly
... ... @@ -541,9 +577,9 @@ void glut_passive_motion(int x, int y) {
541 577  
542 578 if (simulation || build_inlet_outlet && !mods) {
543 579 bool flag = flow.epsilon_edge((float)posX, (float)posY, (float)posZ, eps, direction_index);
544   - if (flag && !flow_direction)
  580 + if (flag && !glyph_mode)
545 581 render_direction = true;
546   - else if (!flag && !flow_direction) {
  582 + else if (!flag && !glyph_mode) {
547 583 if (render_direction) // if the direction is displaying currently, do a short delay
548 584 Sleep(300);
549 585 render_direction = false;
... ... @@ -654,6 +690,7 @@ void glut_mouse(int button, int state, int x, int y) {
654 690 fragment_ratio = 0.0f;
655 691  
656 692 change_fragment = false;
  693 + render_flow_rate = true;
657 694 flow.modify_synthetic_connection(u, rou, hilbert_curve, height_threshold, in, out, fragment_ratio, default_radius);
658 695 }
659 696 // move connections along y-axis
... ... @@ -799,7 +836,7 @@ void glut_wheel(int wheel, int direction, int x, int y) {
799 836  
800 837 if (!to_select_pressure) {
801 838 bool flag = flow.epsilon_vertex((float)posX, (float)posY, (float)posZ, eps, scale, pressure_index);
802   - if (flag && simulation) {
  839 + if (flag && simulation && !glyph_mode) {
803 840 float tmp_r;
804 841 if (direction > 0) { // increase radii
805 842 tmp_r = flow.get_radius(pressure_index);
... ... @@ -812,6 +849,7 @@ void glut_wheel(int wheel, int direction, int x, int y) {
812 849 tmp_r = default_radius;
813 850 }
814 851 flow.set_radius(pressure_index, tmp_r);
  852 + undo = true; // undo rendering
815 853 }
816 854 else if (!mods) {
817 855 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) {
839 877 else
840 878 scale = 1.0f;
841 879 }
  880 + undo = true;
  881 + redisplay = true;
842 882 }
843 883  
844 884 glutPostRedisplay();
... ... @@ -855,22 +895,44 @@ void glut_keyboard(unsigned char key, int x, int y) {
855 895 flow.save_network();
856 896 break;
857 897  
  898 + // convert network to binary format (.nwt)
  899 + case 'c': {
  900 + std::vector<std::string> tmp = stim::parser::split(args.arg(0), '.');
  901 + std::stringstream ss;
  902 + ss << tmp[0] << ".nwt";
  903 + std::string filename = ss.str();
  904 + flow.saveNwt(filename);
  905 + break;
  906 + }
  907 +
858 908 // flow vector field visualization, Glyphs
859 909 case 'f':
860   - if (flow_direction && !manufacture && (simulation || build_inlet_outlet)) {
861   - flow_direction = false;
  910 + if (glyph_mode && !manufacture && (simulation || build_inlet_outlet)) {
  911 + glyph_mode = false;
  912 + frame_mode = false;
862 913 redisplay = true; // lines and arrows rendering use the same display list
863 914 int num = glutGet(GLUT_MENU_NUM_ITEMS);
864 915 if (num == 1)
865 916 glut_set_menu(num, 2);
866 917 }
867   - else if (!flow_direction && !manufacture && (simulation || build_inlet_outlet)) {
868   - flow_direction = true;
  918 + else if (!glyph_mode && !manufacture && (simulation || build_inlet_outlet)) {
  919 + glyph_mode = true;
869 920 redisplay = true;
870 921 int num = glutGet(GLUT_MENU_NUM_ITEMS);
871 922 if (num == 2)
872 923 glut_set_menu(num, 1);
873 924 }
  925 + undo = true;
  926 + break;
  927 +
  928 + // filaments around arrows
  929 + case 'g':
  930 + if (glyph_mode) {
  931 + if (frame_mode)
  932 + frame_mode = false;
  933 + else
  934 + frame_mode = true;
  935 + }
874 936 break;
875 937  
876 938 // open/close index marks
... ... @@ -879,6 +941,7 @@ void glut_keyboard(unsigned char key, int x, int y) {
879 941 mark_index = false;
880 942 else
881 943 mark_index = true;
  944 + undo = true;
882 945 break;
883 946  
884 947 // output image stack
... ... @@ -947,8 +1010,10 @@ void advertise() {
947 1010  
948 1011 std::cout << "Usage(keyboard): e -> open/close indexing" << std::endl;
949 1012 std::cout << " m -> build synthetic connections(connection mode)/output augmented network as image stack (manufacture mode)" << std::endl;
950   - std::cout << " s -> save network flow profiles in result folder as cvs files" << std::endl;
951   - std::cout << " f -> open/close vector field visualization" << std::endl;
  1013 + std::cout << " s -> save network flow profiles in profile folder as cvs files" << std::endl;
  1014 + std::cout << " c -> convert .obj network to .nwt network and stores in main folder" << std::endl;
  1015 + std::cout << " f -> open/close vector field visualization mode" << std::endl;
  1016 + std::cout << " g -> render filament frames in vector fiedl visualization mode" << std::endl;
952 1017  
953 1018 std::cout << args.str();
954 1019 }
... ... @@ -966,6 +1031,8 @@ int main(int argc, char* argv[]) {
966 1031 args.add("stack", "load the image stack");
967 1032 args.add("stackres", "spacing between pixel samples in each dimension(in units/pixel)", "1 1 1", "real value > 0");
968 1033 args.add("stackdir", "set the directory of the output image stack", "", "any existing directory (ex. /home/name/network)");
  1034 + args.add("scale", "scale down rendering fibers");
  1035 + args.add("lcc", "extract the largest connected component");
969 1036  
970 1037 args.parse(argc, argv); // parse the command line
971 1038  
... ... @@ -983,6 +1050,8 @@ int main(int argc, char* argv[]) {
983 1050 std::vector<std::string> tmp = stim::parser::split(args.arg(0), '.');
984 1051 if ("obj" == tmp[1])
985 1052 flow.load_obj(args.arg(0));
  1053 + else if ("nwt" == tmp[1]) // stim network binary format
  1054 + flow.loadNwt(args.arg(0));
986 1055 else if ("swc" == tmp[1])
987 1056 flow.load_swc(args.arg(0));
988 1057 else {
... ... @@ -990,7 +1059,8 @@ int main(int argc, char* argv[]) {
990 1059 std::exit(1);
991 1060 }
992 1061 }
993   - get_background();
  1062 +
  1063 + // extract the largest connected component
994 1064  
995 1065 // get the units to work on
996 1066 units = args["units"].as_string();
... ... @@ -1016,22 +1086,19 @@ int main(int argc, char* argv[]) {
1016 1086 S.load_images(args["stack"].as_string());
1017 1087 // binary transformation
1018 1088 #ifdef __CUDACC__
1019   - unsigned N = S.samples(); // number of pixels loaded
1020   - unsigned char* d_S; // image stack stored in device
  1089 + unsigned N = S.samples(); // number of pixels loaded
  1090 + unsigned char* d_S; // image stack stored in device
1021 1091 unsigned char* h_S = (unsigned char*)malloc(N * sizeof(unsigned char)); // image stack stored in host
1022 1092 cudaMalloc((void**)&d_S, N * sizeof(unsigned char));
1023 1093 cudaMemcpy(d_S, S.data(), N * sizeof(unsigned char), cudaMemcpyHostToDevice);
1024 1094  
1025   -
1026 1095 unsigned thread = 1024;
1027 1096 unsigned block = N / thread + 1;
1028   - binary_transform <<<block, thread>>> (N, d_S, binary_threshold);
  1097 + binary_transform <<<block, thread>>> (N, d_S, binary_threshold); // binaryzation
1029 1098  
1030 1099 cudaMemcpy(h_S, d_S, N * sizeof(unsigned char), cudaMemcpyDeviceToHost);
1031 1100  
1032 1101 S.copy(h_S);
1033   -
1034   - S.save_images("image????.bmp");
1035 1102 #endif
1036 1103 }
1037 1104  
... ... @@ -1044,6 +1111,10 @@ int main(int argc, char* argv[]) {
1044 1111 if (args["stackdir"].is_set())
1045 1112 stackdir = args["stackdir"].as_string();
1046 1113  
  1114 + // get the scale-down factor is provided
  1115 + if (args["scale"].is_set())
  1116 + scale = args["scale"].as_float();
  1117 +
1047 1118 // glut main loop
1048 1119 bb = flow.boundingbox();
1049 1120 glut_initialize();
... ...