codebase / stimlib

Commit 1599ae65702af1249455c2e98c09ec5f247d1752

Authored by Jiaming Guo
1 parent 4f63a044

add swc.h and functions for loading networks from swc files

Showing 4 changed files with 331 additions and 6 deletions   Show diff stats
Inline Side-by-side
stim/biomodels/centerline.h
  Show/Hide comments   View file @1599ae6
... ... @@ -194,8 +194,8 @@ public:
194 194  
195 195 ////resample a fiber in the network
196 196 stim::centerline<T> resample(T spacing)
197   - {
198   - std::cout<<"fiber::resample()"<<std::endl;
  197 + {
  198 + //std::cout<<"fiber::resample()"<<std::endl;
199 199  
200 200 stim::vec3<T> v; //v-direction vector of the segment
201 201 stim::vec3<T> p; //- intermediate point to be added
... ...
stim/biomodels/network.h
  Show/Hide comments   View file @1599ae6
... ... @@ -10,6 +10,7 @@
10 10 #include <math.h>
11 11 #include <stim/math/vec3.h>
12 12 #include <stim/visualization/obj.h>
  13 +#include <stim/visualization/swc.h>
13 14 #include <stim/visualization/cylinder.h>
14 15 #include <stim/structures/kdtree.cuh>
15 16 #include <stim/cuda/cudatools/timer.h>
... ... @@ -142,7 +143,8 @@ class network{
142 143 protected:
143 144  
144 145 std::vector<edge> E; //list of edges
145   - std::vector<vertex> V; //list of vertices.
  146 + std::vector<vertex> V; //list of vertices.
  147 + std::vector<int> NT; //stores a list of neuronal type for each point in the centerline (will set value only when centerline is built from swc file)
146 148  
147 149 public:
148 150  
... ... @@ -388,6 +390,65 @@ public:
388 390 }
389 391 }
390 392  
  393 + void load_swc(std::string filename) {
  394 + stim::swc<T> S; // create swc variable
  395 + S.load(filename); // load the node information
  396 + S.create_tree(); // link those node according to their linking relationships as a tree
  397 +
  398 + NT.push_back(S.node[0].type); // set the neuronal_type value to the first vertex in the network
  399 + std::vector<unsigned> id2vert; // this list stores the SWC vertex ID associated with each network vertex
  400 + unsigned i[2]; // temporary, IDs associated with the first and last points
  401 + for (unsigned int l = 1; l < S.numL(); l++) { // for every node starts from second one
  402 + NT.push_back(S.node[l].type);
  403 + stim::centerline<T> c3(2); // every fiber contains two vertices
  404 + int id = S.node[l].parent_idx - 1;
  405 + c3[0] = S.node[id].point; // store the begin vertex
  406 + c3[1] = S.node[l].point; // store the end vertex
  407 +
  408 + c3.update(); // update the L information
  409 +
  410 + edge new_edge(c3); // contains two points
  411 +
  412 + //get the first and last vertex IDs for the line
  413 + i[0] = S.node[id].idx; //get the SWC ID for the start point
  414 + i[1] = S.node[l].idx; //get the SWC ID for the end point
  415 +
  416 + std::vector<unsigned>::iterator it; //create an iterator for searching the id2vert array
  417 + unsigned it_idx; //create an integer for the id2vert entry index
  418 +
  419 + //find out if the nodes for this fiber have already been created
  420 + it = find(id2vert.begin(), id2vert.end(), i[0]); //look for the first node
  421 + if (it == id2vert.end()) { //if i[0] hasn't already been used
  422 + vertex new_vertex = new_edge[0]; //create a new vertex, assign it a position
  423 + new_vertex.e[0].push_back(E.size()); //add the current edge as outgoing
  424 + new_edge.v[0] = V.size(); //add the new edge to the edge
  425 + V.push_back(new_vertex); //add the new vertex to the vertex list
  426 + id2vert.push_back(i[0]); //add the ID to the ID->vertex conversion list
  427 + }
  428 + else { //if the vertex already exists
  429 + it_idx = std::distance(id2vert.begin(), it);
  430 + V[it_idx].e[0].push_back(E.size()); //add the current edge as outgoing
  431 + new_edge.v[0] = it_idx;
  432 + }
  433 +
  434 + it = find(id2vert.begin(), id2vert.end(), i[1]); //look for the second ID
  435 + if (it == id2vert.end()) { //if i[1] hasn't already been used
  436 + vertex new_vertex = new_edge[1]; //create a new vertex, assign it a position
  437 + new_vertex.e[1].push_back(E.size()); //add the current edge as incoming
  438 + new_edge.v[1] = V.size(); //add the new vertex to the edge
  439 + V.push_back(new_vertex); //add the new vertex to the vertex list
  440 + id2vert.push_back(i[1]); //add the ID to the ID->vertex conversion list
  441 + }
  442 + else { //if the vertex already exists
  443 + it_idx = std::distance(id2vert.begin(), it);
  444 + V[it_idx].e[1].push_back(E.size()); //add the current edge as incoming
  445 + new_edge.v[1] = it_idx;
  446 + }
  447 +
  448 + E.push_back(new_edge); //push the edge to the list
  449 + }
  450 + }
  451 +
391 452 /// Output the network as a string
392 453 std::string str(){
393 454  
... ...
stim/visualization/gl_network.h
  Show/Hide comments   View file @1599ae6
... ... @@ -107,6 +107,88 @@ public:
107 107 glCallList(dlist); // render the display list
108 108 }
109 109  
  110 + /// render the network centerline from swc file as a series of strips in different colors based on the neuronal type
  111 + /// glCenterline0_swc is for only one input
  112 + void glCenterline0_swc() {
  113 + if (!glIsList(dlist)) { // if dlist isn't a display list, create it
  114 + dlist = glGenLists(1); // generate a display list
  115 + glNewList(dlist, GL_COMPILE); // start a new display list
  116 + for (unsigned e = 0; e < E.size(); e++) {
  117 + int type = NT[e];
  118 + switch (type) {
  119 + case 0:
  120 + glColor3f(1.0f, 1.0f, 1.0f); // white for undefined
  121 + glBegin(GL_LINE_STRIP);
  122 + for (unsigned p = 0; p < E[e].size(); p++) {
  123 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  124 + }
  125 + glEnd();
  126 + break;
  127 + case 1:
  128 + glColor3f(1.0f, 0.0f, 0.0f); // red for soma
  129 + glBegin(GL_LINE_STRIP);
  130 + for (unsigned p = 0; p < E[e].size(); p++) {
  131 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  132 + }
  133 + glEnd();
  134 + break;
  135 + case 2:
  136 + glColor3f(1.0f, 0.5f, 0.0f); // orange for axon
  137 + glBegin(GL_LINE_STRIP);
  138 + for (unsigned p = 0; p < E[e].size(); p++) {
  139 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  140 + }
  141 + glEnd();
  142 + break;
  143 + case 3:
  144 + glColor3f(1.0f, 1.0f, 0.0f); // yellow for undefined
  145 + glBegin(GL_LINE_STRIP);
  146 + for (unsigned p = 0; p < E[e].size(); p++) {
  147 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  148 + }
  149 + glEnd();
  150 + break;
  151 + case 4:
  152 + glColor3f(0.0f, 1.0f, 0.0f); // green for undefined
  153 + glBegin(GL_LINE_STRIP);
  154 + for (unsigned p = 0; p < E[e].size(); p++) {
  155 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  156 + }
  157 + glEnd();
  158 + break;
  159 + case 5:
  160 + glColor3f(0.0f, 1.0f, 1.0f); // verdant for undefined
  161 + glBegin(GL_LINE_STRIP);
  162 + for (unsigned p = 0; p < E[e].size(); p++) {
  163 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  164 + }
  165 + glEnd();
  166 + break;
  167 + case 6:
  168 + glColor3f(0.0f, 0.0f, 1.0f); // blue for undefined
  169 + glBegin(GL_LINE_STRIP);
  170 + for (unsigned p = 0; p < E[e].size(); p++) {
  171 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  172 + }
  173 + glEnd();
  174 + break;
  175 + case 7:
  176 + glColor3f(0.5f, 0.0f, 1.0f); // purple for undefined
  177 + glBegin(GL_LINE_STRIP);
  178 + for (unsigned p = 0; p < E[e].size(); p++) {
  179 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  180 + }
  181 + glEnd();
  182 + break;
  183 + }
  184 + }
  185 + glEndList(); //end the display list
  186 + }
  187 + glCallList(dlist); // render the display list
  188 + }
  189 +
  190 + ///render the network centerline as a series of line strips
  191 + ///colors is based on metric values
110 192 void glCenterline(){
111 193  
112 194 if(!glIsList(dlist)){ //if dlist isn't a display list, create it
... ... @@ -117,7 +199,7 @@ public:
117 199 glBegin(GL_LINE_STRIP);
118 200 for(unsigned p = 0; p < E[e].size(); p++){ //for each point on that edge
119 201 glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); //set the vertex position based on the current point
120   - glTexCoord1f(E[e].r(p)); //set the texture coordinate based on the specified magnitude index
  202 + glTexCoord1f(E[e].r(p)); //set the texture coordinate based on the specified magnitude index
121 203 }
122 204 glEnd();
123 205 }
... ... @@ -140,7 +222,7 @@ public:
140 222 for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge
141 223 stim::circle<T> C1 = E[e].circ(p - 1);
142 224 stim::circle<T> C2 = E[e].circ(p);
143   - C1.set_R(10); // scale the circle to the same
  225 + C1.set_R(10); // scale the circle to the same
144 226 C2.set_R(10);
145 227 std::vector< stim::vec3<T> >Cp1 = C1.points(20);
146 228 std::vector< stim::vec3<T> >Cp2 = C2.points(20);
... ... @@ -149,7 +231,7 @@ public:
149 231 }
150 232 else {
151 233 glColor3f(1.0f, 1.0f, 1.0f); // white color for the un-mapping edges
152   - for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge
  234 + for (unsigned p = 1; p < E[e].size(); p++) { // for each point on that edge
153 235 stim::circle<T> C1 = E[e].circ(p - 1);
154 236 stim::circle<T> C2 = E[e].circ(p);
155 237 C1.set_R(10); // scale the circle to the same
... ...
stim/visualization/swc.h 0 → 100644
  Show/Hide comments   View file @1599ae6
  1 +#ifndef STIM_SWC_H
  2 +#define STIM_SWC_H
  3 +
  4 +#include <vector>
  5 +#include <fstream>
  6 +#include <sstream>
  7 +#include <iostream>
  8 +
  9 +//STIM includes
  10 +#include <stim/math/vec3.h>
  11 +#include <stim/parser/parser.h>
  12 +
  13 +namespace stim {
  14 + namespace swc_tree {
  15 + template <typename T>
  16 + class swc_node {
  17 +
  18 + protected:
  19 + enum neuronal_type { SWC_UNDEFINED, SWC_SOMA, SWC_AXON, SWC_DENDRITE, SWC_APICAL_DENDRITE, SWC_FORK_POINT, SWC_END_POINT, SWC_CUSTOM }; // eight types
  20 + enum node_information { INTEGER_LABEL, NEURONAL_TYPE, X_COORDINATE, Y_COORDINATE, Z_COORDINATE, RADIUS, PARENT_LABEL };
  21 +
  22 + public:
  23 + int idx; // the index of current node start from 1(ambiguity)
  24 + neuronal_type type; // the type of neuronal segmemt
  25 + stim::vec3<T> point; // the point coordinates
  26 + T radius; // the radius at current node
  27 + int parent_idx; // parent idx of current node, -1 when it is origin(soma)
  28 + int level; // tree level
  29 + std::vector<int> son_idx; // son idx of current node
  30 +
  31 + swc_node() { // default constructor
  32 + idx = -1; // set to default -1
  33 + parent_idx = -1; // set to origin -1
  34 + level = -1; // set to default -1
  35 + type = SWC_UNDEFINED; // set to undefined type
  36 + radius = 0; // set to 0
  37 + }
  38 +
  39 + void get_node(std::string line) { // get information from the node point we got
  40 +
  41 + // create a vector to store the information of one node point
  42 + std::vector<std::string> p = stim::parser::split(line, ' ');
  43 +
  44 + // for each information contained in the node point we got
  45 + for (unsigned int i = 0; i < p.size(); i++) {
  46 + std::stringstream ss(p[i]); // create a stringstream object for casting
  47 +
  48 + // store different information
  49 + switch (i) {
  50 + case INTEGER_LABEL:
  51 + ss >> idx; // cast the stringstream to the correct numerical value
  52 + break;
  53 + case NEURONAL_TYPE:
  54 + int tmp_type;
  55 + ss >> tmp_type; // cast the stringstream to the correct numerical value
  56 + type = (neuronal_type)tmp_type;
  57 + break;
  58 + case X_COORDINATE:
  59 + T tmp_X;
  60 + ss >> tmp_X; // cast the stringstream to the correct numerical value
  61 + point[0] = tmp_X; // store the X_coordinate in vec3[0]
  62 + break;
  63 + case Y_COORDINATE:
  64 + T tmp_Y;
  65 + ss >> tmp_Y; // cast the stringstream to the correct numerical value
  66 + point[1] = tmp_Y; // store the Y_coordinate in vec3[1]
  67 + break;
  68 + case Z_COORDINATE:
  69 + T tmp_Z;
  70 + ss >> tmp_Z; // cast the stringstream to the correct numerical value
  71 + point[2] = tmp_Z; // store the Z_coordinate in vec3[2]
  72 + break;
  73 + case RADIUS:
  74 + ss >> radius; // cast the stringstream to the correct numerical value
  75 + break;
  76 + case PARENT_LABEL:
  77 + ss >> parent_idx; // cast the stringstream to the correct numerical value
  78 + break;
  79 + }
  80 + }
  81 + }
  82 + };
  83 + } // end of namespace swc_tree
  84 +
  85 + template <typename T>
  86 + class swc {
  87 + public:
  88 + std::vector< typename swc_tree::swc_node<T>> node;
  89 + swc() {}; // default constructor
  90 +
  91 + // load the swc as tree nodes
  92 + void load(std::string filename) {
  93 +
  94 + // load the file
  95 + std::ifstream infile(filename.c_str());
  96 +
  97 + // if the file is invalid, throw an error
  98 + if (!infile) {
  99 + std::cerr << "STUN::SWC Error loading file" << filename << std::endl;
  100 + exit(-1);
  101 + }
  102 +
  103 + std::string line;
  104 + // skip comment
  105 + while (getline(infile, line)) {
  106 + if ('#' == line[0]) // if it is comment line
  107 + continue;
  108 + else
  109 + break;
  110 + }
  111 +
  112 + unsigned int l = 0; // number of nodes
  113 +
  114 + // get rid of the first/origin node
  115 + swc_tree::swc_node<T> new_node;
  116 + new_node.get_node(line);
  117 + l++;
  118 + node.push_back(new_node); // push back the first node
  119 +
  120 + getline(infile, line); // get a new line
  121 + // keep reading the following node point information as string
  122 + while (!line.empty()) { // test for the last empty line
  123 + l++; // still remaining node to be read
  124 +
  125 + swc_tree::swc_node<T> next_node;
  126 + next_node.get_node(line);
  127 + node.push_back(next_node);
  128 +
  129 + getline(infile, line); // get a new line
  130 + }
  131 + }
  132 +
  133 + // read the head comment from swc file
  134 + void read_comment(std::string filename) {
  135 +
  136 + // load the file
  137 + std::ifstream infile(filename.c_str());
  138 +
  139 + // if the file is invalid, throw an error
  140 + if (!infile) {
  141 + std::cerr << "STUN::SWC Error loading file" << filename << std::endl;
  142 + exit(1);
  143 + }
  144 +
  145 + std::string line;
  146 + while (getline(infile, line)) {
  147 + if ('#' == line[0]) {
  148 + std::cout << line << std::endl; // print the comment line by line
  149 + }
  150 + else
  151 + break; // break when reaches to node information
  152 + }
  153 + }
  154 +
  155 + // link those nodes to create a tree
  156 + void create_tree() {
  157 + unsigned n = node.size(); // get the total number of node point
  158 + int cur_level = 0;
  159 +
  160 + // build the origin(soma)
  161 + node[0].level = cur_level;
  162 +
  163 + // go through follow nodes
  164 + for (unsigned i = 1; i < n; i++) {
  165 + if (node[i].parent_idx != node[i - 1].parent_idx)
  166 + cur_level = node[node[i].parent_idx - 1].level + 1;
  167 + node[i].level = cur_level;
  168 + int tmp_parent_idx = node[i].parent_idx - 1; // get the parent node loop idx of current node
  169 + node[tmp_parent_idx].son_idx.push_back(i + 1); // son_idx stores the real idx = loop idx + 1
  170 + }
  171 + }
  172 +
  173 + // return the number of point in swc
  174 + unsigned int numL() {
  175 + return node.size();
  176 + }
  177 +
  178 +
  179 + };
  180 +} // end of namespace stim
  181 +
  182 +#endif
0 183 \ No newline at end of file
... ...