From 586740361645287847389e57d555cd0bc208bd6b Mon Sep 17 00:00:00 2001 From: Jiaming Guo Date: Thu, 12 Jan 2017 11:45:23 -0600 Subject: [PATCH] fix minor error for loading and rendering networks from swc files --- stim/biomodels/network.h | 22 +++++++++++++++------- stim/visualization/cylinder.h | 33 +++++++++++++++++++++------------ stim/visualization/gl_network.h | 179 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------- stim/visualization/swc.h | 12 ++++++------ 4 files changed, 192 insertions(+), 54 deletions(-) diff --git a/stim/biomodels/network.h b/stim/biomodels/network.h index 438bac5..d1b33c0 100644 --- a/stim/biomodels/network.h +++ b/stim/biomodels/network.h @@ -390,6 +390,7 @@ public: } } + //load a network from an SWC file void load_swc(std::string filename) { stim::swc S; // create swc variable S.load(filename); // load the node information @@ -401,16 +402,23 @@ public: for (unsigned int l = 1; l < S.numL(); l++) { // for every node starts from second one NT.push_back(S.node[l].type); stim::centerline c3(2); // every fiber contains two vertices - int id = S.node[l].parent_idx - 1; - c3[0] = S.node[id].point; // store the begin vertex + int p_idx = S.node[l].parent_idx - 1; // get the parent node loop idx of current node + c3[0] = S.node[p_idx].point; // store the begin vertex c3[1] = S.node[l].point; // store the end vertex c3.update(); // update the L information - - edge new_edge(c3); // contains two points - + + stim::cylinder C3(c3); // create a new cylinder in order to copy the origin radius information + // upadate the R information + std::vector radius; + radius.push_back(S.node[p_idx].radius); + radius.push_back(S.node[l].radius); + C3.copy_r(radius); + + edge new_edge(C3); // contains two points + //get the first and last vertex IDs for the line - i[0] = S.node[id].idx; //get the SWC ID for the start point + i[0] = S.node[p_idx].idx; //get the SWC ID for the start point i[1] = S.node[l].idx; //get the SWC ID for the end point std::vector::iterator it; //create an iterator for searching the id2vert array @@ -471,7 +479,7 @@ public: stim::network resample(T spacing){ stim::network n; //create a new network that will be an exact copy, with resampled fibers n.V = V; //copy all vertices - + n.NT = NT; //copy all the neuronal type information n.E.resize(edges()); //allocate space for the edge list //copy all fibers, resampling them in the process diff --git a/stim/visualization/cylinder.h b/stim/visualization/cylinder.h index 4635a34..12b21ac 100644 --- a/stim/visualization/cylinder.h +++ b/stim/visualization/cylinder.h @@ -54,6 +54,12 @@ public: //add_mag(r); } + //copy the original radius + void copy_r(std::vector radius) { + for (unsigned i = 0; i < radius.size(); i++) + R[i] = radius[i]; + } + ///Returns magnitude i at the given length into the fiber (based on the pvalue). ///Interpolates the position along the line. ///@param l: the location of the in the cylinder. @@ -175,20 +181,23 @@ public: /// Integrates a magnitude value along the cylinder. /// @param m is the magnitude value to be integrated (this is usually the radius) T integrate() { + T sum = 0; //initialize the integral to zero + if (L.size() == 1) + return sum; + else { + T m0, m1; //allocate space for both magnitudes in a single segment + m0 = R[0]; //initialize the first point and magnitude to the first point in the cylinder + T len = L[1]; //allocate space for the segment length - T sum = 0; //initialize the integral to zero - T m0, m1; //allocate space for both magnitudes in a single segment - m0 = R[0]; //initialize the first point and magnitude to the first point in the cylinder - T len = L[1]; //allocate space for the segment length - - - for (unsigned p = 1; p < size(); p++) { //for every consecutive point in the cylinder - m1 = R[p]; - if (p > 1) len = (L[p] - L[p - 1]); //calculate the segment length using the L array - sum += (m0 + m1) / (T)2.0 * len; //add the average magnitude, weighted by the segment length - m0 = m1; //move to the next segment by shifting points + + for (unsigned p = 1; p < size(); p++) { //for every consecutive point in the cylinder + m1 = R[p]; + if (p > 1) len = (L[p] - L[p - 1]); //calculate the segment length using the L array + sum += (m0 + m1) / (T)2.0 * len; //add the average magnitude, weighted by the segment length + m0 = m1; //move to the next segment by shifting points + } + return sum; //return the integral } - return sum; //return the integral } /// Resamples the cylinder to provide a maximum distance of "spacing" between centerline points. All current diff --git a/stim/visualization/gl_network.h b/stim/visualization/gl_network.h index e1028c0..0a7f172 100644 --- a/stim/visualization/gl_network.h +++ b/stim/visualization/gl_network.h @@ -43,28 +43,6 @@ public: return bb; //return the bounding box } - //void renderCylinder(T x1, T y1, T z1, T x2, T y2, T z2, T radius, int subdivisions) { - // T dx = x2 - x1; - // T dy = y2 - y1; - // T dz = z2 - z1; - // /// handle the degenerate case with an approximation - // if (dz == 0) - // dz = .00000001; - // T d = sqrt(dx*dx + dy*dy + dz*dz); - // T ax = 57.2957795*acos(dz / d); // 180°/pi - // if (dz < 0.0) - // ax = -ax; - // T rx = -dy*dz; - // T ry = dx*dz; - - // glPushMatrix(); - // glTranslatef(x1, y1, z1); - // glRotatef(ax, rx, ry, 0.0); - - // glutSolidCylinder(radius, d, subdivisions, 1); - // glPopMatrix(); - //} - ///render cylinder based on points from the top/bottom hat ///@param C1 set of points from one of the hat void renderCylinder(std::vector< stim::vec3 > C1, std::vector< stim::vec3 > C2) { @@ -114,7 +92,7 @@ public: dlist = glGenLists(1); // generate a display list glNewList(dlist, GL_COMPILE); // start a new display list for (unsigned e = 0; e < E.size(); e++) { - int type = NT[e]; + int type = NT[e]; // get the neuronal type switch (type) { case 0: glColor3f(1.0f, 1.0f, 1.0f); // white for undefined @@ -187,8 +165,8 @@ public: glCallList(dlist); // render the display list } - ///render the network centerline as a series of line strips - ///colors is based on metric values + ///render the network centerline as a series of line strips(when loading at least two networks, otherwise using glCenterline0()) + ///colors are based on metric values void glCenterline(){ if(!glIsList(dlist)){ //if dlist isn't a display list, create it @@ -208,11 +186,39 @@ public: glCallList(dlist); //render the display list } + ///render the network cylinder as a series of tubes + ///colors are based on metric values + void glCylinder() { + if (!glIsList(dlist)) { //if dlist isn't a display list, create it + dlist = glGenLists(1); //generate a display list + glNewList(dlist, GL_COMPILE); //start a new display list + for (unsigned e = 0; e < E.size(); e++) { //for each edge in the network + for (unsigned p = 1; p < E[e].size(); p++) { // for each point on that edge + stim::circle C1 = E[e].circ(p - 1); + stim::circle C2 = E[e].circ(p); + C1.set_R(10); // scale the circle to the same + C2.set_R(10); + std::vector< stim::vec3 >Cp1 = C1.points(20); + std::vector< stim::vec3 >Cp2 = C2.points(20); + glBegin(GL_QUAD_STRIP); + for (unsigned i = 0; i < Cp1.size(); i++) { // for every point on the circle + glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]); + glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]); + glTexCoord1f(E[e].r(p)); + } + glEnd(); + } //set the texture coordinate based on the specified magnitude index + } + glEndList(); //end the display list + } + glCallList(dlist); //render the display list + } + ///render the GT network cylinder as series of tubes ///@param dlist1 is the display list ///@param map is the mapping relationship between two networks ///@param colormap is the random generated color set for render - void glCylinderGT(GLuint &dlist1, std::vector map, std::vector colormap) { + void glRandColorCylinder1(GLuint &dlist1, std::vector map, std::vector colormap) { if (!glIsList(dlist1)) { // if dlist1 isn't a display list, create it dlist1 = glGenLists(1); // generate a display list glNewList(dlist1, GL_COMPILE); // start a new display list @@ -258,11 +264,57 @@ public: glCallList(dlist1); } + void glRandColorCylinder1_swc(GLuint &dlist1, std::vector map, std::vector colormap) { + if (!glIsList(dlist1)) { // if dlist1 isn't a display list, create it + dlist1 = glGenLists(1); // generate a display list + glNewList(dlist1, GL_COMPILE); // start a new display list + for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network + if (map[e] != unsigned(-1)) { + glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]); + for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge + stim::circle C1 = E[e].circ(p - 1); + stim::circle C2 = E[e].circ(p); + C1.set_R(0.5); // scale the circle to the same + C2.set_R(0.5); + std::vector< stim::vec3 >Cp1 = C1.points(20); + std::vector< stim::vec3 >Cp2 = C2.points(20); + renderCylinder(Cp1, Cp2); + } + } + else { + glColor3f(1.0f, 1.0f, 1.0f); // white color for the un-mapping edges + for (unsigned p = 1; p < E[e].size(); p++) { // for each point on that edge + stim::circle C1 = E[e].circ(p - 1); + stim::circle C2 = E[e].circ(p); + C1.set_R(0.5); // scale the circle to the same + C2.set_R(0.5); + std::vector< stim::vec3 >Cp1 = C1.points(20); + std::vector< stim::vec3 >Cp2 = C2.points(20); + renderCylinder(Cp1, Cp2); + } + } + } + for (unsigned v = 0; v < V.size(); v++) { + size_t num_edge = V[v].e[0].size() + V[v].e[1].size(); + if (num_edge > 1) { // if it is the joint vertex + glColor3f(0.3, 0.3, 0.3); // gray color + renderBall(V[v][0], V[v][1], V[v][2], 1, 20); + } + else { // if it is the terminal vertex + glColor3f(0.6, 0.6, 0.6); // more white gray + renderBall(V[v][0], V[v][1], V[v][2], 1, 20); + } + } + glEndList(); + } + glCallList(dlist1); + } + ///render the T network cylinder as series of tubes ///@param dlist2 is the display list ///@param map is the mapping relationship between two networks ///@param colormap is the random generated color set for render - void glCylinderT(GLuint &dlist2, std::vector map, std::vector colormap) { + void glRandColorCylinder2(GLuint &dlist2, std::vector map, std::vector colormap) { if (!glIsList(dlist2)) { dlist2 = glGenLists(1); glNewList(dlist2, GL_COMPILE); @@ -308,11 +360,57 @@ public: glCallList(dlist2); } + void glRandColorCylinder2_swc(GLuint &dlist2, std::vector map, std::vector colormap) { + if (!glIsList(dlist2)) { + dlist2 = glGenLists(1); + glNewList(dlist2, GL_COMPILE); + for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network + if (map[e] != unsigned(-1)) { + glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]); + for (unsigned p = 0; p < E[e].size() - 1; p++) {// for each point on that edge + stim::circle C1 = E[e].circ(p); + stim::circle C2 = E[e].circ(p + 1); + C1.set_R(0.5); // scale the circle to the same + C2.set_R(0.5); + std::vector< stim::vec3 >Cp1 = C1.points(20); + std::vector< stim::vec3 >Cp2 = C2.points(20); + renderCylinder(Cp1, Cp2); + } + } + else { + glColor3f(1.0f, 1.0f, 1.0f); // white color for the un-mapping edges + for (unsigned p = 0; p < E[e].size() - 1; p++) {// for each point on that edge + stim::circle C1 = E[e].circ(p); + stim::circle C2 = E[e].circ(p + 1); + C1.set_R(0.5); // scale the circle to the same + C2.set_R(0.5); + std::vector< stim::vec3 >Cp1 = C1.points(20); + std::vector< stim::vec3 >Cp2 = C2.points(20); + renderCylinder(Cp1, Cp2); + } + } + } + for (unsigned v = 0; v < V.size(); v++) { + size_t num_edge = V[v].e[0].size() + V[v].e[1].size(); + if (num_edge > 1) { // if it is the joint vertex + glColor3f(0.3, 0.3, 0.3); // gray color + renderBall(V[v][0], V[v][1], V[v][2], 1, 20); + } + else { // if it is the terminal vertex + glColor3f(0.6, 0.6, 0.6); // more white gray + renderBall(V[v][0], V[v][1], V[v][2], 1, 20); + } + } + glEndList(); + } + glCallList(dlist2); + } + /// Render the GT network centerline as a series of line strips in random different color ///@param dlist1 is the display list ///@param map is the mapping relationship between two networks ///@param colormap is the random generated color set for render - void glRandColorCenterlineGT(GLuint &dlist1, std::vector map, std::vector colormap) { + void glRandColorCenterline1(GLuint &dlist1, std::vector map, std::vector colormap) { if (!glIsList(dlist1)) { dlist1 = glGenLists(1); glNewList(dlist1, GL_COMPILE); @@ -343,7 +441,7 @@ public: ///@param dlist2 is the display list ///@param map is the mapping relationship between two networks ///@param colormap is the random generated color set for render - void glRandColorCenterlineT(GLuint &dlist2, std::vector map, std::vector colormap) { + void glRandColorCenterline2(GLuint &dlist2, std::vector map, std::vector colormap) { if (!glIsList(dlist2)) { dlist2 = glGenLists(1); glNewList(dlist2, GL_COMPILE); @@ -444,6 +542,29 @@ public: // glCallList(dlist2); //} + + //void renderCylinder(T x1, T y1, T z1, T x2, T y2, T z2, T radius, int subdivisions) { + // T dx = x2 - x1; + // T dy = y2 - y1; + // T dz = z2 - z1; + // /// handle the degenerate case with an approximation + // if (dz == 0) + // dz = .00000001; + // T d = sqrt(dx*dx + dy*dy + dz*dz); + // T ax = 57.2957795*acos(dz / d); // 180°/pi + // if (dz < 0.0) + // ax = -ax; + // T rx = -dy*dz; + // T ry = dx*dz; + + // glPushMatrix(); + // glTranslatef(x1, y1, z1); + // glRotatef(ax, rx, ry, 0.0); + + // glutSolidCylinder(radius, d, subdivisions, 1); + // glPopMatrix(); + //} + }; //end stim::gl_network class }; //end stim namespace diff --git a/stim/visualization/swc.h b/stim/visualization/swc.h index 8166f38..beaffe0 100644 --- a/stim/visualization/swc.h +++ b/stim/visualization/swc.h @@ -20,7 +20,7 @@ namespace stim { enum node_information { INTEGER_LABEL, NEURONAL_TYPE, X_COORDINATE, Y_COORDINATE, Z_COORDINATE, RADIUS, PARENT_LABEL }; public: - int idx; // the index of current node start from 1(ambiguity) + int idx; // the index of current node start from 1(original index from the file) neuronal_type type; // the type of neuronal segmemt stim::vec3 point; // the point coordinates T radius; // the radius at current node @@ -44,7 +44,7 @@ namespace stim { // for each information contained in the node point we got for (unsigned int i = 0; i < p.size(); i++) { std::stringstream ss(p[i]); // create a stringstream object for casting - + // store different information switch (i) { case INTEGER_LABEL: @@ -103,8 +103,8 @@ namespace stim { std::string line; // skip comment while (getline(infile, line)) { - if ('#' == line[0]) // if it is comment line - continue; + if ('#' == line[0] || line.empty()) // if it is comment line or empty line + continue; // keep read else break; } @@ -144,7 +144,7 @@ namespace stim { std::string line; while (getline(infile, line)) { - if ('#' == line[0]) { + if ('#' == line[0] || line.empty()) { std::cout << line << std::endl; // print the comment line by line } else @@ -165,7 +165,7 @@ namespace stim { if (node[i].parent_idx != node[i - 1].parent_idx) cur_level = node[node[i].parent_idx - 1].level + 1; node[i].level = cur_level; - int tmp_parent_idx = node[i].parent_idx - 1; // get the parent node loop idx of current node + int tmp_parent_idx = node[i].parent_idx - 1; // get the parent node loop idx of current node node[tmp_parent_idx].son_idx.push_back(i + 1); // son_idx stores the real idx = loop idx + 1 } } -- libgit2 0.21.4