diff --git a/stim/biomodels/centerline.h b/stim/biomodels/centerline.h index 7a4ca45..136aab7 100644 --- a/stim/biomodels/centerline.h +++ b/stim/biomodels/centerline.h @@ -194,8 +194,8 @@ public: ////resample a fiber in the network stim::centerline resample(T spacing) - { - std::cout<<"fiber::resample()"< v; //v-direction vector of the segment stim::vec3 p; //- intermediate point to be added diff --git a/stim/biomodels/network.h b/stim/biomodels/network.h index 3354a4d..438bac5 100644 --- a/stim/biomodels/network.h +++ b/stim/biomodels/network.h @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -142,7 +143,8 @@ class network{ protected: std::vector E; //list of edges - std::vector V; //list of vertices. + std::vector V; //list of vertices. + std::vector NT; //stores a list of neuronal type for each point in the centerline (will set value only when centerline is built from swc file) public: @@ -388,6 +390,65 @@ public: } } + void load_swc(std::string filename) { + 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 + + 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 + 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 + c3[1] = S.node[l].point; // store the end vertex + + c3.update(); // update the L information + + 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[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 + unsigned it_idx; //create an integer for the id2vert entry index + + //find out if the nodes for this fiber have already been created + it = find(id2vert.begin(), id2vert.end(), i[0]); //look for the first node + if (it == id2vert.end()) { //if i[0] hasn't already been used + vertex new_vertex = new_edge[0]; //create a new vertex, assign it a position + new_vertex.e[0].push_back(E.size()); //add the current edge as outgoing + new_edge.v[0] = V.size(); //add the new edge to the edge + V.push_back(new_vertex); //add the new vertex to the vertex list + id2vert.push_back(i[0]); //add the ID to the ID->vertex conversion list + } + else { //if the vertex already exists + it_idx = std::distance(id2vert.begin(), it); + V[it_idx].e[0].push_back(E.size()); //add the current edge as outgoing + new_edge.v[0] = it_idx; + } + + 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 + 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 + id2vert.push_back(i[1]); //add the ID to the ID->vertex conversion list + } + else { //if the vertex already exists + it_idx = std::distance(id2vert.begin(), it); + V[it_idx].e[1].push_back(E.size()); //add the current edge as incoming + new_edge.v[1] = it_idx; + } + + E.push_back(new_edge); //push the edge to the list + } + } + /// Output the network as a string std::string str(){ diff --git a/stim/visualization/gl_network.h b/stim/visualization/gl_network.h index ce4c00a..e1028c0 100644 --- a/stim/visualization/gl_network.h +++ b/stim/visualization/gl_network.h @@ -107,6 +107,88 @@ public: glCallList(dlist); // render the display list } + /// render the network centerline from swc file as a series of strips in different colors based on the neuronal type + /// glCenterline0_swc is for only one input + void glCenterline0_swc() { + 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++) { + int type = NT[e]; + switch (type) { + case 0: + glColor3f(1.0f, 1.0f, 1.0f); // white for undefined + 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(); + break; + case 1: + glColor3f(1.0f, 0.0f, 0.0f); // red for soma + 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(); + break; + case 2: + glColor3f(1.0f, 0.5f, 0.0f); // orange for axon + 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(); + break; + case 3: + glColor3f(1.0f, 1.0f, 0.0f); // yellow for undefined + 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(); + break; + case 4: + glColor3f(0.0f, 1.0f, 0.0f); // green for undefined + 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(); + break; + case 5: + glColor3f(0.0f, 1.0f, 1.0f); // verdant for undefined + 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(); + break; + case 6: + glColor3f(0.0f, 0.0f, 1.0f); // blue for undefined + 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(); + break; + case 7: + glColor3f(0.5f, 0.0f, 1.0f); // purple for undefined + 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(); + break; + } + } + glEndList(); //end the display list + } + glCallList(dlist); // render the display list + } + + ///render the network centerline as a series of line strips + ///colors is based on metric values void glCenterline(){ if(!glIsList(dlist)){ //if dlist isn't a display list, create it @@ -117,7 +199,7 @@ public: glBegin(GL_LINE_STRIP); for(unsigned p = 0; p < E[e].size(); p++){ //for each point on that edge glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); //set the vertex position based on the current point - glTexCoord1f(E[e].r(p)); //set the texture coordinate based on the specified magnitude index + glTexCoord1f(E[e].r(p)); //set the texture coordinate based on the specified magnitude index } glEnd(); } @@ -140,7 +222,7 @@ public: 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 + 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); @@ -149,7 +231,7 @@ public: } 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 + 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 diff --git a/stim/visualization/swc.h b/stim/visualization/swc.h new file mode 100644 index 0000000..8166f38 --- /dev/null +++ b/stim/visualization/swc.h @@ -0,0 +1,182 @@ +#ifndef STIM_SWC_H +#define STIM_SWC_H + +#include +#include +#include +#include + +//STIM includes +#include +#include + +namespace stim { + namespace swc_tree { + template + class swc_node { + + protected: + enum neuronal_type { SWC_UNDEFINED, SWC_SOMA, SWC_AXON, SWC_DENDRITE, SWC_APICAL_DENDRITE, SWC_FORK_POINT, SWC_END_POINT, SWC_CUSTOM }; // eight types + 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) + neuronal_type type; // the type of neuronal segmemt + stim::vec3 point; // the point coordinates + T radius; // the radius at current node + int parent_idx; // parent idx of current node, -1 when it is origin(soma) + int level; // tree level + std::vector son_idx; // son idx of current node + + swc_node() { // default constructor + idx = -1; // set to default -1 + parent_idx = -1; // set to origin -1 + level = -1; // set to default -1 + type = SWC_UNDEFINED; // set to undefined type + radius = 0; // set to 0 + } + + void get_node(std::string line) { // get information from the node point we got + + // create a vector to store the information of one node point + std::vector p = stim::parser::split(line, ' '); + + // 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: + ss >> idx; // cast the stringstream to the correct numerical value + break; + case NEURONAL_TYPE: + int tmp_type; + ss >> tmp_type; // cast the stringstream to the correct numerical value + type = (neuronal_type)tmp_type; + break; + case X_COORDINATE: + T tmp_X; + ss >> tmp_X; // cast the stringstream to the correct numerical value + point[0] = tmp_X; // store the X_coordinate in vec3[0] + break; + case Y_COORDINATE: + T tmp_Y; + ss >> tmp_Y; // cast the stringstream to the correct numerical value + point[1] = tmp_Y; // store the Y_coordinate in vec3[1] + break; + case Z_COORDINATE: + T tmp_Z; + ss >> tmp_Z; // cast the stringstream to the correct numerical value + point[2] = tmp_Z; // store the Z_coordinate in vec3[2] + break; + case RADIUS: + ss >> radius; // cast the stringstream to the correct numerical value + break; + case PARENT_LABEL: + ss >> parent_idx; // cast the stringstream to the correct numerical value + break; + } + } + } + }; + } // end of namespace swc_tree + + template + class swc { + public: + std::vector< typename swc_tree::swc_node> node; + swc() {}; // default constructor + + // load the swc as tree nodes + void load(std::string filename) { + + // load the file + std::ifstream infile(filename.c_str()); + + // if the file is invalid, throw an error + if (!infile) { + std::cerr << "STUN::SWC Error loading file" << filename << std::endl; + exit(-1); + } + + std::string line; + // skip comment + while (getline(infile, line)) { + if ('#' == line[0]) // if it is comment line + continue; + else + break; + } + + unsigned int l = 0; // number of nodes + + // get rid of the first/origin node + swc_tree::swc_node new_node; + new_node.get_node(line); + l++; + node.push_back(new_node); // push back the first node + + getline(infile, line); // get a new line + // keep reading the following node point information as string + while (!line.empty()) { // test for the last empty line + l++; // still remaining node to be read + + swc_tree::swc_node next_node; + next_node.get_node(line); + node.push_back(next_node); + + getline(infile, line); // get a new line + } + } + + // read the head comment from swc file + void read_comment(std::string filename) { + + // load the file + std::ifstream infile(filename.c_str()); + + // if the file is invalid, throw an error + if (!infile) { + std::cerr << "STUN::SWC Error loading file" << filename << std::endl; + exit(1); + } + + std::string line; + while (getline(infile, line)) { + if ('#' == line[0]) { + std::cout << line << std::endl; // print the comment line by line + } + else + break; // break when reaches to node information + } + } + + // link those nodes to create a tree + void create_tree() { + unsigned n = node.size(); // get the total number of node point + int cur_level = 0; + + // build the origin(soma) + node[0].level = cur_level; + + // go through follow nodes + for (unsigned i = 1; i < n; i++) { + 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 + node[tmp_parent_idx].son_idx.push_back(i + 1); // son_idx stores the real idx = loop idx + 1 + } + } + + // return the number of point in swc + unsigned int numL() { + return node.size(); + } + + + }; +} // end of namespace stim + +#endif \ No newline at end of file -- libgit2 0.21.4