diff --git a/stim/biomodels/network.h b/stim/biomodels/network.h index 782f418..8b804b4 100644 --- a/stim/biomodels/network.h +++ b/stim/biomodels/network.h @@ -395,31 +395,40 @@ public: stim::swc S; // create swc variable S.load(filename); // load the node information S.create_tree(); // link those node according to their linking relationships as a tree + S.resample(); //NT.push_back(S.node[0].type); // set the neuronal_type value to the first vertex in the network std::vector id2vert; // this list stores the SWC vertex ID associated with each network vertex unsigned i[2]; // temporary, IDs associated with the first and last points - for (unsigned int l = 1; l < S.numL(); l++) { // for every node starts from second one + + for (unsigned int l = 0; l < S.numE(); l++) { // for every edge //NT.push_back(S.node[l].type); - stim::centerline c3(2); // every fiber contains two vertices - 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 + + std::vector< stim::vec3 > c; + S.get_points(l, c); + + stim::centerline c3(c.size()); // new fiber + + for (unsigned j = 0; j < c.size(); j++) + c3[j] = c[j]; // copy the points c3.update(); // update the L information 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); + S.get_radius(l, radius); + C3.copy_r(radius); - edge new_edge(C3); // contains two points + edge new_edge(C3); // new edge + + //create an edge from the given centerline + unsigned int I = new_edge.size(); //calculate the number of points on the centerline //get the first and last vertex IDs for the line - 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 + i[0] = S.E[l].front(); + i[1] = S.E[l].back(); std::vector::iterator it; //create an iterator for searching the id2vert array unsigned it_idx; //create an integer for the id2vert entry index @@ -441,7 +450,7 @@ public: it = find(id2vert.begin(), id2vert.end(), i[1]); //look for the second ID if (it == id2vert.end()) { //if i[1] hasn't already been used - vertex new_vertex = new_edge[1]; //create a new vertex, assign it a position + vertex new_vertex = new_edge[I - 1]; //create a new vertex, assign it a position new_vertex.e[1].push_back(E.size()); //add the current edge as incoming new_edge.v[1] = V.size(); //add the new vertex to the edge V.push_back(new_vertex); //add the new vertex to the vertex list @@ -620,8 +629,6 @@ public: kdt.create(c, n_data, MaxTreeLevels); for(unsigned e = 0; e < R.E.size(); e++){ //for each edge in A - R.E[e].add_mag(0); //add a new magnitude for the metric - size_t errormag_id = R.E[e].nmags() - 1; size_t n = R.E[e].size(); // the number of points in current edge T* query = new T[3 * n]; @@ -634,8 +641,8 @@ public: kdt.cpu_search(queryPt, n, nnIdx, dists); //find the distance between A and the current network for(unsigned p = 0; p < R.E[e].size(); p++){ - m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma); //calculate the metric value based on the distance - R.E[e].set_mag(errormag_id, p, m1[p]); //set the error for the second point in the segment + m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma); //calculate the metric value based on the distance + R.E[e].set_r(p, m1[p]); //set the error for the second point in the segment } } #endif @@ -746,10 +753,8 @@ public: edge_point_num[ee] = B.E[ee].size(); for(unsigned e = 0; e < A.E.size(); e++){ //for each edge in A - A.E[e].add_mag(0); //add a new magnitude for the metric - size_t errormag_id = A.E[e].nmags() - 1; - - size_t n = A.E[e].size(); // the number of edges in A + + size_t n = A.E[e].size(); //the number of edges in A T* queryPt = new T[3 * n]; T* m1 = new T[n]; @@ -761,7 +766,7 @@ public: for(unsigned p = 0; p < A.E[e].size(); p++){ m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma); //calculate the metric value based on the distance - A.E[e].set_mag(errormag_id, p, m1[p]); //set the error for the second point in the segment + A.E[e].set_r(p, m1[p]); //set the error for the second point in the segment unsigned id = 0; //mapping edge's idx size_t num = 0; //total number of points before #th edge @@ -909,9 +914,8 @@ public: for(unsigned e = 0; e < R.E.size(); e++){ // for each edge in A T M; // the sum of metrics of current edge - unsigned errormag_id = R.E[e].nmags() - 1; for(unsigned p = 0; p < R.E[e].size(); p++) - M += A.E[e].m(p, errormag_id); + M += A.E[e].r(p); M = M / A.E[e].size(); if(M > threshold) C[e] = (unsigned)-1; diff --git a/stim/visualization/gl_network.h b/stim/visualization/gl_network.h index e6ef107..3107256 100644 --- a/stim/visualization/gl_network.h +++ b/stim/visualization/gl_network.h @@ -115,6 +115,8 @@ public: 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 + glEnable(GL_COLOR_MATERIAL); + glColorMaterial(GL_FRONT, GL_AMBIENT_AND_DIFFUSE); 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); @@ -125,7 +127,7 @@ public: std::vector< stim::vec3 >Cp2 = C2.glpoints(20); glBegin(GL_QUAD_STRIP); for (unsigned i = 0; i < Cp1.size(); i++) { // for every point on the circle(+1 means closing the circle) - glTexCoord1f(E[e].r(p-1)); + glTexCoord1f(E[e].r(p - 1)); glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]); glTexCoord1f(E[e].r(p)); glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]); @@ -150,77 +152,67 @@ public: 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 glRandColorCylinder1(GLuint &dlist1, std::vector map, std::vector colormap, float sigma, float radius) { + ///render a T as a adjoint network of GT in transparancy(darkgreen, overlaid) + void glAdjointCylinder(float sigma, float radius) { - if (radius != sigma) // if render radius was changed by user, create a new display list - glDeleteLists(dlist1, 1); - 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(2*radius); // scale the circle to the same - C2.set_R(2*radius); - std::vector< stim::vec3 >Cp1 = C1.glpoints(20); - std::vector< stim::vec3 >Cp2 = C2.glpoints(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(2*radius); // scale the circle to the same - C2.set_R(2*radius); - std::vector< stim::vec3 >Cp1 = C1.glpoints(20); - std::vector< stim::vec3 >Cp2 = C2.glpoints(20); - renderCylinder(Cp1, Cp2); + if (radius != sigma) // if render radius was changed by user, create a new display list + glDeleteLists(dlist+4, 1); + if (!glIsList(dlist+4)) { + glNewList(dlist+4, GL_COMPILE); + 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(2 * radius); // scale the circle to the same + C2.set_R(2 * radius); + std::vector< stim::vec3 >Cp1 = C1.glpoints(20); + std::vector< stim::vec3 >Cp2 = C2.glpoints(20); + glBegin(GL_QUAD_STRIP); + for (unsigned i = 0; i < Cp1.size(); i++) { // for every point on the circle(+1 means closing the circle) + glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]); + glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]); } - } + glEnd(); + } // set the texture coordinate based on the specified magnitude index } - 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], 3*radius, 20); + for (unsigned n = 0; n < V.size(); n++) { + size_t num = V[n].e[0].size(); // store the number of outgoing edge of that vertex + if (num != 0) { // if it has outgoing edge + unsigned idx = V[n].e[0][0]; // find the index of first outgoing edge of that vertex } - 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], 3*radius, 20); + else { + unsigned idx = V[n].e[1][0]; // find the index of first incoming edge of that vertex } + renderBall(V[n][0], V[n][1], V[n][2], 2 * radius, 20); } glEndList(); } - glCallList(dlist1); + glCallList(dlist+4); } - ///render the T network cylinder as series of tubes - ///@param dlist2 is the display list + ///render the network cylinder as series of tubes + ///@param I is a indicator: 0 -> GT, 1 -> T ///@param map is the mapping relationship between two networks ///@param colormap is the random generated color set for render - void glRandColorCylinder2(GLuint &dlist2, std::vector map, std::vector colormap, float sigma, float radius) { + void glRandColorCylinder(int I, std::vector map, std::vector colormap, float sigma, float radius) { - if (radius != sigma) // if render radius was changed by user, create a new display list - glDeleteLists(dlist2, 1); - 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 (radius != sigma) // if render radius was changed by user, create a new display list + glDeleteLists(dlist+2, 1); + if (!glIsList(dlist+2)) { // if dlist isn't a display list, create it + glNewList(dlist+2, GL_COMPILE); // start a new display list + glEnable(GL_COLOR_MATERIAL); + glColorMaterial(GL_FRONT, GL_AMBIENT_AND_DIFFUSE); + 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 = 1; p < E[e].size(); p++) { // for each point on that edge(except last one) + if(I == 0) // if it is to render GT + glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]); + else // if it is to render T + glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[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(2*radius); // scale the circle to the same + C1.set_R(2*radius); // scale the circle to the same C2.set_R(2*radius); std::vector< stim::vec3 >Cp1 = C1.glpoints(20); std::vector< stim::vec3 >Cp2 = C2.glpoints(20); @@ -252,52 +244,27 @@ public: } } glEndList(); + glDisable(GL_COLOR_MATERIAL); } - glCallList(dlist2); + glCallList(dlist+2); } - /// Render the GT network centerline as a series of line strips in random different color + ///render the network centerline as a series of line strips in random different color + ///@param I is a indicator: 0 -> GT, 1 -> T ///@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 glRandColorCenterline1(GLuint &dlist1, std::vector map, std::vector colormap) { + void glRandColorCenterline(int I, GLuint &dlist1, std::vector map, std::vector colormap) { if (!glIsList(dlist1)) { dlist1 = glGenLists(1); glNewList(dlist1, GL_COMPILE); for (unsigned e = 0; e < E.size(); e++) { - if (map[e] != unsigned(-1)) { // if it has corresponding edge in another network - glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]); - glBegin(GL_LINE_STRIP); - for (unsigned p = 0; p < E[e].size(); p++) { - glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); - } - glEnd(); - } - else { - glColor3f(1.0, 1.0, 1.0); // white color - glBegin(GL_LINE_STRIP); - for (unsigned p = 0; p < E[e].size(); p++) { - glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); - } - glEnd(); - } - } - glEndList(); - } - glCallList(dlist1); - } + if (map[e] != unsigned(-1)) { // if it has corresponding edge in another network + if (I == 0) // if it is to render GT + glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]); + else // if it is to render T + glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]); - /// Render the T network centerline as a series of line strips in random different color - ///@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 glRandColorCenterline2(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++) { - if (map[e] != unsigned(-1)) { // if it has corresponding edge in another network - glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]); glBegin(GL_LINE_STRIP); for (unsigned p = 0; p < E[e].size(); p++) { glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); @@ -305,7 +272,7 @@ public: glEnd(); } else { - glColor3f(1.0, 1.0, 1.0); // white color + glColor3f(1.0, 1.0, 1.0); // white color glBegin(GL_LINE_STRIP); for (unsigned p = 0; p < E[e].size(); p++) { glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); @@ -315,7 +282,7 @@ public: } glEndList(); } - glCallList(dlist2); + glCallList(dlist1); } //void glRandColorCenterlineGT(GLuint &dlist1, std::vector map, std::vector colormap){ diff --git a/stim/visualization/swc.h b/stim/visualization/swc.h index 634ab5b..8a61edb 100644 --- a/stim/visualization/swc.h +++ b/stim/visualization/swc.h @@ -86,6 +86,8 @@ namespace stim { class swc { public: std::vector< typename swc_tree::swc_node > node; + std::vector< std::vector > E; // list of nodes + swc() {}; // default constructor // load the swc as tree nodes @@ -170,11 +172,104 @@ namespace stim { } } + // create a new edge + void create_up(std::vector& edge,swc_tree::swc_node cur_node, int target) { + while (cur_node.parent_idx != target) { + edge.push_back(cur_node.idx); + cur_node = node[cur_node.parent_idx - 1]; // move to next node + } + edge.push_back(cur_node.idx); // push back the start/end vertex of current edge + + std::reverse(edge.begin(), edge.end()); // follow the original flow direction + } + void create_down(std::vector& edge, swc_tree::swc_node cur_node, int j) { + while (cur_node.son_idx.size() != 0) { + edge.push_back(cur_node.idx); + if (cur_node.son_idx.size() > 1) + cur_node = node[cur_node.son_idx[j] - 1]; // move to next node + else + cur_node = node[cur_node.son_idx[0] - 1]; + } + edge.push_back(cur_node.idx); // push back the start/end vertex of current edge + } + + // resample the tree-like SWC + void resample() { + unsigned n = node.size(); + + std::vector joint_node; + for (unsigned i = 1; i < n; i++) { // search all nodes(except the first one) to find joint nodes + if (node[i].son_idx.size() > 1) { + joint_node.push_back(node[i].idx); // store the original index + } + } + + std::vector new_edge; // new edge in the network + + n = joint_node.size(); + + for (unsigned i = 0; i < n; i++) { // for every joint nodes + std::vector new_edge; // new edge in the network + swc_tree::swc_node cur_node = node[joint_node[i] - 1]; // store current node + + // go up + swc_tree::swc_node tmp_node = node[cur_node.parent_idx - 1]; + while (tmp_node.parent_idx != -1 && tmp_node.son_idx.size() == 1) { + tmp_node = node[tmp_node.parent_idx - 1]; + } + int target = tmp_node.parent_idx; // store the go-up target + create_up(new_edge, cur_node, target); + E.push_back(new_edge); + new_edge.clear(); + + // go down + unsigned t = cur_node.son_idx.size(); + for (unsigned j = 0; j < t; j++) { // for every son node + tmp_node = node[cur_node.son_idx[j] - 1]; // move down + while (tmp_node.son_idx.size() == 1) { + tmp_node = node[tmp_node.son_idx[0] - 1]; + } + if (tmp_node.son_idx.size() == 0) { + create_down(new_edge, cur_node, j); + E.push_back(new_edge); + new_edge.clear(); + } + } + } + } + + // get points in one edge + void get_points(int e, std::vector< stim::vec3 >& V) { + V.resize(E[e].size()); + for (unsigned i = 0; i < E[e].size(); i++) { + unsigned id = E[e][i] - 1; + + V[i][0] = node[id].point[0]; + V[i][1] = node[id].point[1]; + V[i][2] = node[id].point[2]; + } + } + + // get radius information in one edge + void get_radius(int e, std::vector& radius) { + radius.resize(E[e].size()); + for (unsigned i = 0; i < E[e].size(); i++) { + unsigned id = E[e][i] - 1; + + radius[i] = node[id].radius; + } + } + // return the number of point in swc - unsigned int numL() { + unsigned int numP() { return node.size(); } + // return the number of edges in swc after resample + unsigned int numE() { + return E.size(); + } + }; } // end of namespace stim -- libgit2 0.21.4