codebase / stimlib

Commit 586740361645287847389e57d555cd0bc208bd6b

Authored by Jiaming Guo
1 parent 8705dfb0

fix minor error for loading and rendering networks from swc files

Showing 4 changed files with 192 additions and 54 deletions   Show diff stats
Inline Side-by-side
stim/biomodels/network.h
  Show/Hide comments   View file @5867403
... ... @@ -390,6 +390,7 @@ public:
390 390 }
391 391 }
392 392  
  393 + //load a network from an SWC file
393 394 void load_swc(std::string filename) {
394 395 stim::swc<T> S; // create swc variable
395 396 S.load(filename); // load the node information
... ... @@ -401,16 +402,23 @@ public:
401 402 for (unsigned int l = 1; l < S.numL(); l++) { // for every node starts from second one
402 403 NT.push_back(S.node[l].type);
403 404 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
  405 + int p_idx = S.node[l].parent_idx - 1; // get the parent node loop idx of current node
  406 + c3[0] = S.node[p_idx].point; // store the begin vertex
406 407 c3[1] = S.node[l].point; // store the end vertex
407 408  
408 409 c3.update(); // update the L information
409   -
410   - edge new_edge(c3); // contains two points
411   -
  410 +
  411 + stim::cylinder<T> C3(c3); // create a new cylinder in order to copy the origin radius information
  412 + // upadate the R information
  413 + std::vector<T> radius;
  414 + radius.push_back(S.node[p_idx].radius);
  415 + radius.push_back(S.node[l].radius);
  416 + C3.copy_r(radius);
  417 +
  418 + edge new_edge(C3); // contains two points
  419 +
412 420 //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
  421 + i[0] = S.node[p_idx].idx; //get the SWC ID for the start point
414 422 i[1] = S.node[l].idx; //get the SWC ID for the end point
415 423  
416 424 std::vector<unsigned>::iterator it; //create an iterator for searching the id2vert array
... ... @@ -471,7 +479,7 @@ public:
471 479 stim::network<T> resample(T spacing){
472 480 stim::network<T> n; //create a new network that will be an exact copy, with resampled fibers
473 481 n.V = V; //copy all vertices
474   -
  482 + n.NT = NT; //copy all the neuronal type information
475 483 n.E.resize(edges()); //allocate space for the edge list
476 484  
477 485 //copy all fibers, resampling them in the process
... ...
stim/visualization/cylinder.h
  Show/Hide comments   View file @5867403
... ... @@ -54,6 +54,12 @@ public:
54 54 //add_mag(r);
55 55 }
56 56  
  57 + //copy the original radius
  58 + void copy_r(std::vector<T> radius) {
  59 + for (unsigned i = 0; i < radius.size(); i++)
  60 + R[i] = radius[i];
  61 + }
  62 +
57 63 ///Returns magnitude i at the given length into the fiber (based on the pvalue).
58 64 ///Interpolates the position along the line.
59 65 ///@param l: the location of the in the cylinder.
... ... @@ -175,20 +181,23 @@ public:
175 181 /// Integrates a magnitude value along the cylinder.
176 182 /// @param m is the magnitude value to be integrated (this is usually the radius)
177 183 T integrate() {
  184 + T sum = 0; //initialize the integral to zero
  185 + if (L.size() == 1)
  186 + return sum;
  187 + else {
  188 + T m0, m1; //allocate space for both magnitudes in a single segment
  189 + m0 = R[0]; //initialize the first point and magnitude to the first point in the cylinder
  190 + T len = L[1]; //allocate space for the segment length
178 191  
179   - T sum = 0; //initialize the integral to zero
180   - T m0, m1; //allocate space for both magnitudes in a single segment
181   - m0 = R[0]; //initialize the first point and magnitude to the first point in the cylinder
182   - T len = L[1]; //allocate space for the segment length
183   -
184   -
185   - for (unsigned p = 1; p < size(); p++) { //for every consecutive point in the cylinder
186   - m1 = R[p];
187   - if (p > 1) len = (L[p] - L[p - 1]); //calculate the segment length using the L array
188   - sum += (m0 + m1) / (T)2.0 * len; //add the average magnitude, weighted by the segment length
189   - m0 = m1; //move to the next segment by shifting points
  192 +
  193 + for (unsigned p = 1; p < size(); p++) { //for every consecutive point in the cylinder
  194 + m1 = R[p];
  195 + if (p > 1) len = (L[p] - L[p - 1]); //calculate the segment length using the L array
  196 + sum += (m0 + m1) / (T)2.0 * len; //add the average magnitude, weighted by the segment length
  197 + m0 = m1; //move to the next segment by shifting points
  198 + }
  199 + return sum; //return the integral
190 200 }
191   - return sum; //return the integral
192 201 }
193 202  
194 203 /// Resamples the cylinder to provide a maximum distance of "spacing" between centerline points. All current
... ...
stim/visualization/gl_network.h
  Show/Hide comments   View file @5867403
... ... @@ -43,28 +43,6 @@ public:
43 43 return bb; //return the bounding box
44 44 }
45 45  
46   - //void renderCylinder(T x1, T y1, T z1, T x2, T y2, T z2, T radius, int subdivisions) {
47   - // T dx = x2 - x1;
48   - // T dy = y2 - y1;
49   - // T dz = z2 - z1;
50   - // /// handle the degenerate case with an approximation
51   - // if (dz == 0)
52   - // dz = .00000001;
53   - // T d = sqrt(dx*dx + dy*dy + dz*dz);
54   - // T ax = 57.2957795*acos(dz / d); // 180°/pi
55   - // if (dz < 0.0)
56   - // ax = -ax;
57   - // T rx = -dy*dz;
58   - // T ry = dx*dz;
59   -
60   - // glPushMatrix();
61   - // glTranslatef(x1, y1, z1);
62   - // glRotatef(ax, rx, ry, 0.0);
63   -
64   - // glutSolidCylinder(radius, d, subdivisions, 1);
65   - // glPopMatrix();
66   - //}
67   -
68 46 ///render cylinder based on points from the top/bottom hat
69 47 ///@param C1 set of points from one of the hat
70 48 void renderCylinder(std::vector< stim::vec3<T> > C1, std::vector< stim::vec3<T> > C2) {
... ... @@ -114,7 +92,7 @@ public:
114 92 dlist = glGenLists(1); // generate a display list
115 93 glNewList(dlist, GL_COMPILE); // start a new display list
116 94 for (unsigned e = 0; e < E.size(); e++) {
117   - int type = NT[e];
  95 + int type = NT[e]; // get the neuronal type
118 96 switch (type) {
119 97 case 0:
120 98 glColor3f(1.0f, 1.0f, 1.0f); // white for undefined
... ... @@ -187,8 +165,8 @@ public:
187 165 glCallList(dlist); // render the display list
188 166 }
189 167  
190   - ///render the network centerline as a series of line strips
191   - ///colors is based on metric values
  168 + ///render the network centerline as a series of line strips(when loading at least two networks, otherwise using glCenterline0())
  169 + ///colors are based on metric values
192 170 void glCenterline(){
193 171  
194 172 if(!glIsList(dlist)){ //if dlist isn't a display list, create it
... ... @@ -208,11 +186,39 @@ public:
208 186 glCallList(dlist); //render the display list
209 187 }
210 188  
  189 + ///render the network cylinder as a series of tubes
  190 + ///colors are based on metric values
  191 + void glCylinder() {
  192 + if (!glIsList(dlist)) { //if dlist isn't a display list, create it
  193 + dlist = glGenLists(1); //generate a display list
  194 + glNewList(dlist, GL_COMPILE); //start a new display list
  195 + for (unsigned e = 0; e < E.size(); e++) { //for each edge in the network
  196 + for (unsigned p = 1; p < E[e].size(); p++) { // for each point on that edge
  197 + stim::circle<T> C1 = E[e].circ(p - 1);
  198 + stim::circle<T> C2 = E[e].circ(p);
  199 + C1.set_R(10); // scale the circle to the same
  200 + C2.set_R(10);
  201 + std::vector< stim::vec3<T> >Cp1 = C1.points(20);
  202 + std::vector< stim::vec3<T> >Cp2 = C2.points(20);
  203 + glBegin(GL_QUAD_STRIP);
  204 + for (unsigned i = 0; i < Cp1.size(); i++) { // for every point on the circle
  205 + glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]);
  206 + glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]);
  207 + glTexCoord1f(E[e].r(p));
  208 + }
  209 + glEnd();
  210 + } //set the texture coordinate based on the specified magnitude index
  211 + }
  212 + glEndList(); //end the display list
  213 + }
  214 + glCallList(dlist); //render the display list
  215 + }
  216 +
211 217 ///render the GT network cylinder as series of tubes
212 218 ///@param dlist1 is the display list
213 219 ///@param map is the mapping relationship between two networks
214 220 ///@param colormap is the random generated color set for render
215   - void glCylinderGT(GLuint &dlist1, std::vector<unsigned> map, std::vector<T> colormap) {
  221 + void glRandColorCylinder1(GLuint &dlist1, std::vector<unsigned> map, std::vector<T> colormap) {
216 222 if (!glIsList(dlist1)) { // if dlist1 isn't a display list, create it
217 223 dlist1 = glGenLists(1); // generate a display list
218 224 glNewList(dlist1, GL_COMPILE); // start a new display list
... ... @@ -258,11 +264,57 @@ public:
258 264 glCallList(dlist1);
259 265 }
260 266  
  267 + void glRandColorCylinder1_swc(GLuint &dlist1, std::vector<unsigned> map, std::vector<T> colormap) {
  268 + if (!glIsList(dlist1)) { // if dlist1 isn't a display list, create it
  269 + dlist1 = glGenLists(1); // generate a display list
  270 + glNewList(dlist1, GL_COMPILE); // start a new display list
  271 + for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network
  272 + if (map[e] != unsigned(-1)) {
  273 + glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]);
  274 + for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge
  275 + stim::circle<T> C1 = E[e].circ(p - 1);
  276 + stim::circle<T> C2 = E[e].circ(p);
  277 + C1.set_R(0.5); // scale the circle to the same
  278 + C2.set_R(0.5);
  279 + std::vector< stim::vec3<T> >Cp1 = C1.points(20);
  280 + std::vector< stim::vec3<T> >Cp2 = C2.points(20);
  281 + renderCylinder(Cp1, Cp2);
  282 + }
  283 + }
  284 + else {
  285 + glColor3f(1.0f, 1.0f, 1.0f); // white color for the un-mapping edges
  286 + for (unsigned p = 1; p < E[e].size(); p++) { // for each point on that edge
  287 + stim::circle<T> C1 = E[e].circ(p - 1);
  288 + stim::circle<T> C2 = E[e].circ(p);
  289 + C1.set_R(0.5); // scale the circle to the same
  290 + C2.set_R(0.5);
  291 + std::vector< stim::vec3<T> >Cp1 = C1.points(20);
  292 + std::vector< stim::vec3<T> >Cp2 = C2.points(20);
  293 + renderCylinder(Cp1, Cp2);
  294 + }
  295 + }
  296 + }
  297 + for (unsigned v = 0; v < V.size(); v++) {
  298 + size_t num_edge = V[v].e[0].size() + V[v].e[1].size();
  299 + if (num_edge > 1) { // if it is the joint vertex
  300 + glColor3f(0.3, 0.3, 0.3); // gray color
  301 + renderBall(V[v][0], V[v][1], V[v][2], 1, 20);
  302 + }
  303 + else { // if it is the terminal vertex
  304 + glColor3f(0.6, 0.6, 0.6); // more white gray
  305 + renderBall(V[v][0], V[v][1], V[v][2], 1, 20);
  306 + }
  307 + }
  308 + glEndList();
  309 + }
  310 + glCallList(dlist1);
  311 + }
  312 +
261 313 ///render the T network cylinder as series of tubes
262 314 ///@param dlist2 is the display list
263 315 ///@param map is the mapping relationship between two networks
264 316 ///@param colormap is the random generated color set for render
265   - void glCylinderT(GLuint &dlist2, std::vector<unsigned> map, std::vector<T> colormap) {
  317 + void glRandColorCylinder2(GLuint &dlist2, std::vector<unsigned> map, std::vector<T> colormap) {
266 318 if (!glIsList(dlist2)) {
267 319 dlist2 = glGenLists(1);
268 320 glNewList(dlist2, GL_COMPILE);
... ... @@ -308,11 +360,57 @@ public:
308 360 glCallList(dlist2);
309 361 }
310 362  
  363 + void glRandColorCylinder2_swc(GLuint &dlist2, std::vector<unsigned> map, std::vector<T> colormap) {
  364 + if (!glIsList(dlist2)) {
  365 + dlist2 = glGenLists(1);
  366 + glNewList(dlist2, GL_COMPILE);
  367 + for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network
  368 + if (map[e] != unsigned(-1)) {
  369 + glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]);
  370 + for (unsigned p = 0; p < E[e].size() - 1; p++) {// for each point on that edge
  371 + stim::circle<T> C1 = E[e].circ(p);
  372 + stim::circle<T> C2 = E[e].circ(p + 1);
  373 + C1.set_R(0.5); // scale the circle to the same
  374 + C2.set_R(0.5);
  375 + std::vector< stim::vec3<T> >Cp1 = C1.points(20);
  376 + std::vector< stim::vec3<T> >Cp2 = C2.points(20);
  377 + renderCylinder(Cp1, Cp2);
  378 + }
  379 + }
  380 + else {
  381 + glColor3f(1.0f, 1.0f, 1.0f); // white color for the un-mapping edges
  382 + for (unsigned p = 0; p < E[e].size() - 1; p++) {// for each point on that edge
  383 + stim::circle<T> C1 = E[e].circ(p);
  384 + stim::circle<T> C2 = E[e].circ(p + 1);
  385 + C1.set_R(0.5); // scale the circle to the same
  386 + C2.set_R(0.5);
  387 + std::vector< stim::vec3<T> >Cp1 = C1.points(20);
  388 + std::vector< stim::vec3<T> >Cp2 = C2.points(20);
  389 + renderCylinder(Cp1, Cp2);
  390 + }
  391 + }
  392 + }
  393 + for (unsigned v = 0; v < V.size(); v++) {
  394 + size_t num_edge = V[v].e[0].size() + V[v].e[1].size();
  395 + if (num_edge > 1) { // if it is the joint vertex
  396 + glColor3f(0.3, 0.3, 0.3); // gray color
  397 + renderBall(V[v][0], V[v][1], V[v][2], 1, 20);
  398 + }
  399 + else { // if it is the terminal vertex
  400 + glColor3f(0.6, 0.6, 0.6); // more white gray
  401 + renderBall(V[v][0], V[v][1], V[v][2], 1, 20);
  402 + }
  403 + }
  404 + glEndList();
  405 + }
  406 + glCallList(dlist2);
  407 + }
  408 +
311 409 /// Render the GT network centerline as a series of line strips in random different color
312 410 ///@param dlist1 is the display list
313 411 ///@param map is the mapping relationship between two networks
314 412 ///@param colormap is the random generated color set for render
315   - void glRandColorCenterlineGT(GLuint &dlist1, std::vector<unsigned> map, std::vector<T> colormap) {
  413 + void glRandColorCenterline1(GLuint &dlist1, std::vector<unsigned> map, std::vector<T> colormap) {
316 414 if (!glIsList(dlist1)) {
317 415 dlist1 = glGenLists(1);
318 416 glNewList(dlist1, GL_COMPILE);
... ... @@ -343,7 +441,7 @@ public:
343 441 ///@param dlist2 is the display list
344 442 ///@param map is the mapping relationship between two networks
345 443 ///@param colormap is the random generated color set for render
346   - void glRandColorCenterlineT(GLuint &dlist2, std::vector<unsigned> map, std::vector<T> colormap) {
  444 + void glRandColorCenterline2(GLuint &dlist2, std::vector<unsigned> map, std::vector<T> colormap) {
347 445 if (!glIsList(dlist2)) {
348 446 dlist2 = glGenLists(1);
349 447 glNewList(dlist2, GL_COMPILE);
... ... @@ -444,6 +542,29 @@ public:
444 542 // glCallList(dlist2);
445 543 //}
446 544  
  545 +
  546 + //void renderCylinder(T x1, T y1, T z1, T x2, T y2, T z2, T radius, int subdivisions) {
  547 + // T dx = x2 - x1;
  548 + // T dy = y2 - y1;
  549 + // T dz = z2 - z1;
  550 + // /// handle the degenerate case with an approximation
  551 + // if (dz == 0)
  552 + // dz = .00000001;
  553 + // T d = sqrt(dx*dx + dy*dy + dz*dz);
  554 + // T ax = 57.2957795*acos(dz / d); // 180°/pi
  555 + // if (dz < 0.0)
  556 + // ax = -ax;
  557 + // T rx = -dy*dz;
  558 + // T ry = dx*dz;
  559 +
  560 + // glPushMatrix();
  561 + // glTranslatef(x1, y1, z1);
  562 + // glRotatef(ax, rx, ry, 0.0);
  563 +
  564 + // glutSolidCylinder(radius, d, subdivisions, 1);
  565 + // glPopMatrix();
  566 + //}
  567 +
447 568 }; //end stim::gl_network class
448 569 }; //end stim namespace
449 570  
... ...
stim/visualization/swc.h
  Show/Hide comments   View file @5867403
... ... @@ -20,7 +20,7 @@ namespace stim {
20 20 enum node_information { INTEGER_LABEL, NEURONAL_TYPE, X_COORDINATE, Y_COORDINATE, Z_COORDINATE, RADIUS, PARENT_LABEL };
21 21  
22 22 public:
23   - int idx; // the index of current node start from 1(ambiguity)
  23 + int idx; // the index of current node start from 1(original index from the file)
24 24 neuronal_type type; // the type of neuronal segmemt
25 25 stim::vec3<T> point; // the point coordinates
26 26 T radius; // the radius at current node
... ... @@ -44,7 +44,7 @@ namespace stim {
44 44 // for each information contained in the node point we got
45 45 for (unsigned int i = 0; i < p.size(); i++) {
46 46 std::stringstream ss(p[i]); // create a stringstream object for casting
47   -
  47 +
48 48 // store different information
49 49 switch (i) {
50 50 case INTEGER_LABEL:
... ... @@ -103,8 +103,8 @@ namespace stim {
103 103 std::string line;
104 104 // skip comment
105 105 while (getline(infile, line)) {
106   - if ('#' == line[0]) // if it is comment line
107   - continue;
  106 + if ('#' == line[0] || line.empty()) // if it is comment line or empty line
  107 + continue; // keep read
108 108 else
109 109 break;
110 110 }
... ... @@ -144,7 +144,7 @@ namespace stim {
144 144  
145 145 std::string line;
146 146 while (getline(infile, line)) {
147   - if ('#' == line[0]) {
  147 + if ('#' == line[0] || line.empty()) {
148 148 std::cout << line << std::endl; // print the comment line by line
149 149 }
150 150 else
... ... @@ -165,7 +165,7 @@ namespace stim {
165 165 if (node[i].parent_idx != node[i - 1].parent_idx)
166 166 cur_level = node[node[i].parent_idx - 1].level + 1;
167 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
  168 + int tmp_parent_idx = node[i].parent_idx - 1; // get the parent node loop idx of current node
169 169 node[tmp_parent_idx].son_idx.push_back(i + 1); // son_idx stores the real idx = loop idx + 1
170 170 }
171 171 }
... ...