codebase / stimlib

Commit 18f207ce4d9af0b33c505a3b5492df63d11f189e

Authored by David Mayerich
2 parents 5bd6c411 2df16996

Merge branch 'JACK' into 'master' 

Jack

See merge request !27
Showing 3 changed files with 233 additions and 120 deletions   Show diff stats
Inline Side-by-side
stim/biomodels/network.h
  Show/Hide comments   View file @18f207c
... ... @@ -395,31 +395,40 @@ public:
395 395 stim::swc<T> S; // create swc variable
396 396 S.load(filename); // load the node information
397 397 S.create_tree(); // link those node according to their linking relationships as a tree
  398 + S.resample();
398 399  
399 400 //NT.push_back(S.node[0].type); // set the neuronal_type value to the first vertex in the network
400 401 std::vector<unsigned> id2vert; // this list stores the SWC vertex ID associated with each network vertex
401 402 unsigned i[2]; // temporary, IDs associated with the first and last points
402   - for (unsigned int l = 1; l < S.numL(); l++) { // for every node starts from second one
  403 +
  404 + for (unsigned int l = 0; l < S.numE(); l++) { // for every edge
403 405 //NT.push_back(S.node[l].type);
404   - stim::centerline<T> c3(2); // every fiber contains two vertices
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
407   - c3[1] = S.node[l].point; // store the end vertex
  406 +
  407 + std::vector< stim::vec3<T> > c;
  408 + S.get_points(l, c);
  409 +
  410 + stim::centerline<T> c3(c.size()); // new fiber
  411 +
  412 + for (unsigned j = 0; j < c.size(); j++)
  413 + c3[j] = c[j]; // copy the points
408 414  
409 415 c3.update(); // update the L information
410 416  
411 417 stim::cylinder<T> C3(c3); // create a new cylinder in order to copy the origin radius information
412 418 // upadate the R information
413 419 std::vector<T> radius;
414   - radius.push_back(S.node[p_idx].radius);
415   - radius.push_back(S.node[l].radius);
  420 + S.get_radius(l, radius);
  421 +
416 422 C3.copy_r(radius);
417 423  
418   - edge new_edge(C3); // contains two points
  424 + edge new_edge(C3); // new edge
  425 +
  426 + //create an edge from the given centerline
  427 + unsigned int I = new_edge.size(); //calculate the number of points on the centerline
419 428  
420 429 //get the first and last vertex IDs for the line
421   - i[0] = S.node[p_idx].idx; //get the SWC ID for the start point
422   - i[1] = S.node[l].idx; //get the SWC ID for the end point
  430 + i[0] = S.E[l].front();
  431 + i[1] = S.E[l].back();
423 432  
424 433 std::vector<unsigned>::iterator it; //create an iterator for searching the id2vert array
425 434 unsigned it_idx; //create an integer for the id2vert entry index
... ... @@ -441,7 +450,7 @@ public:
441 450  
442 451 it = find(id2vert.begin(), id2vert.end(), i[1]); //look for the second ID
443 452 if (it == id2vert.end()) { //if i[1] hasn't already been used
444   - vertex new_vertex = new_edge[1]; //create a new vertex, assign it a position
  453 + vertex new_vertex = new_edge[I - 1]; //create a new vertex, assign it a position
445 454 new_vertex.e[1].push_back(E.size()); //add the current edge as incoming
446 455 new_edge.v[1] = V.size(); //add the new vertex to the edge
447 456 V.push_back(new_vertex); //add the new vertex to the vertex list
... ... @@ -620,8 +629,6 @@ public:
620 629 kdt.create(c, n_data, MaxTreeLevels);
621 630  
622 631 for(unsigned e = 0; e < R.E.size(); e++){ //for each edge in A
623   - R.E[e].add_mag(0); //add a new magnitude for the metric
624   - size_t errormag_id = R.E[e].nmags() - 1;
625 632  
626 633 size_t n = R.E[e].size(); // the number of points in current edge
627 634 T* query = new T[3 * n];
... ... @@ -634,8 +641,8 @@ public:
634 641 kdt.cpu_search(queryPt, n, nnIdx, dists); //find the distance between A and the current network
635 642  
636 643 for(unsigned p = 0; p < R.E[e].size(); p++){
637   - m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma); //calculate the metric value based on the distance
638   - R.E[e].set_mag(errormag_id, p, m1[p]); //set the error for the second point in the segment
  644 + m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma); //calculate the metric value based on the distance
  645 + R.E[e].set_r(p, m1[p]); //set the error for the second point in the segment
639 646 }
640 647 }
641 648 #endif
... ... @@ -746,10 +753,8 @@ public:
746 753 edge_point_num[ee] = B.E[ee].size();
747 754  
748 755 for(unsigned e = 0; e < A.E.size(); e++){ //for each edge in A
749   - A.E[e].add_mag(0); //add a new magnitude for the metric
750   - size_t errormag_id = A.E[e].nmags() - 1;
751   -
752   - size_t n = A.E[e].size(); // the number of edges in A
  756 +
  757 + size_t n = A.E[e].size(); //the number of edges in A
753 758  
754 759 T* queryPt = new T[3 * n];
755 760 T* m1 = new T[n];
... ... @@ -761,7 +766,7 @@ public:
761 766  
762 767 for(unsigned p = 0; p < A.E[e].size(); p++){
763 768 m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma); //calculate the metric value based on the distance
764   - A.E[e].set_mag(errormag_id, p, m1[p]); //set the error for the second point in the segment
  769 + A.E[e].set_r(p, m1[p]); //set the error for the second point in the segment
765 770  
766 771 unsigned id = 0; //mapping edge's idx
767 772 size_t num = 0; //total number of points before #th edge
... ... @@ -909,9 +914,8 @@ public:
909 914  
910 915 for(unsigned e = 0; e < R.E.size(); e++){ // for each edge in A
911 916 T M; // the sum of metrics of current edge
912   - unsigned errormag_id = R.E[e].nmags() - 1;
913 917 for(unsigned p = 0; p < R.E[e].size(); p++)
914   - M += A.E[e].m(p, errormag_id);
  918 + M += A.E[e].r(p);
915 919 M = M / A.E[e].size();
916 920 if(M > threshold)
917 921 C[e] = (unsigned)-1;
... ...
stim/visualization/gl_network.h
  Show/Hide comments   View file @18f207c
... ... @@ -55,7 +55,7 @@ public:
55 55 //glFlush();
56 56 }
57 57  
58   - ///render the vertex as sphere
  58 + ///render the vertex as sphere based on glut build-in function
59 59 ///@param x, y, z are the three coordinates of the center point
60 60 ///@param radius is the radius of the sphere
61 61 ///@param subdivisions is the slice/stride along/around z-axis
... ... @@ -66,6 +66,53 @@ public:
66 66 glPopMatrix();
67 67 }
68 68  
  69 + ///render the vertex as sphere based on transformation
  70 + ///@param x, y, z are the three coordinates of the center point
  71 + ///@param radius is the radius of the sphere
  72 + ///@param slice is the number of subdivisions around the z-axis
  73 + ///@param stack is the number of subdivisions along the z-axis
  74 + void renderBall(T x, T y, T z, T radius, T slice, T stack) {
  75 + T step_z = stim::PI / slice; // step angle along z-axis
  76 + T step_xy = 2 * stim::PI / stack; // step angle in xy-plane
  77 + T xx[4], yy[4], zz[4]; // store coordinates
  78 +
  79 + T angle_z = 0.0; // start angle
  80 + T angle_xy = 0.0;
  81 +
  82 + glBegin(GL_QUADS);
  83 + for (unsigned i = 0; i < slice; i++) // around the z-axis
  84 + {
  85 + angle_z = i * step_z; // step step_z each time
  86 +
  87 + for (unsigned j = 0; j < stack; j++) // along the z-axis
  88 + {
  89 + angle_xy = j * step_xy; // step step_xy each time, draw floor by floor
  90 +
  91 + xx[0] = radius * std::sin(angle_z) * std::cos(angle_xy); // four vertices
  92 + yy[0] = radius * std::sin(angle_z) * std::sin(angle_xy);
  93 + zz[0] = radius * std::cos(angle_z);
  94 +
  95 + xx[1] = radius * std::sin(angle_z + step_z) * std::cos(angle_xy);
  96 + yy[1] = radius * std::sin(angle_z + step_z) * std::sin(angle_xy);
  97 + zz[1] = radius * std::cos(angle_z + step_z);
  98 +
  99 + xx[2] = radius * std::sin(angle_z + step_z) * std::cos(angle_xy + step_xy);
  100 + yy[2] = radius * std::sin(angle_z + step_z) * std::sin(angle_xy + step_xy);
  101 + zz[2] = radius * std::cos(angle_z + step_z);
  102 +
  103 + xx[3] = radius * std::sin(angle_z) * std::cos(angle_xy + step_xy);
  104 + yy[3] = radius * std::sin(angle_z) * std::sin(angle_xy + step_xy);
  105 + zz[3] = radius * std::cos(angle_z);
  106 +
  107 + for (unsigned k = 0; k < 4; k++)
  108 + {
  109 + glVertex3f(x + xx[k], y + yy[k], z + zz[k]); // draw the floor plane
  110 + }
  111 + }
  112 + }
  113 + glEnd();
  114 + }
  115 +
69 116 /// Render the network centerline as a series of line strips.
70 117 /// glCenterline0 is for only one input
71 118 void glCenterline0(){
... ... @@ -115,6 +162,8 @@ public:
115 162 if (!glIsList(dlist)) { // if dlist isn't a display list, create it
116 163 dlist = glGenLists(1); // generate a display list
117 164 glNewList(dlist, GL_COMPILE); // start a new display list
  165 + glEnable(GL_COLOR_MATERIAL);
  166 + glColorMaterial(GL_FRONT, GL_AMBIENT_AND_DIFFUSE);
118 167 for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network
119 168 for (unsigned p = 1; p < E[e].size(); p++) { // for each point on that edge
120 169 stim::circle<T> C1 = E[e].circ(p - 1);
... ... @@ -125,7 +174,7 @@ public:
125 174 std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
126 175 glBegin(GL_QUAD_STRIP);
127 176 for (unsigned i = 0; i < Cp1.size(); i++) { // for every point on the circle(+1 means closing the circle)
128   - glTexCoord1f(E[e].r(p-1));
  177 + glTexCoord1f(E[e].r(p - 1));
129 178 glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]);
130 179 glTexCoord1f(E[e].r(p));
131 180 glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]);
... ... @@ -143,84 +192,74 @@ public:
143 192 unsigned idx = V[n].e[1][0]; // find the index of first incoming edge of that vertex
144 193 glTexCoord1f(E[idx].r(E[idx].size() - 1)); // bind the texture as metric of last point on that edge
145 194 }
146   - renderBall(V[n][0], V[n][1], V[n][2], 2*radius, 20);
  195 + renderBall(V[n][0], V[n][1], V[n][2], 2*radius, 20, 20);
147 196 }
148 197 glEndList(); // end the display list
149 198 }
150 199 glCallList(dlist); // render the display list
151 200 }
152 201  
153   - ///render the GT network cylinder as series of tubes
154   - ///@param dlist1 is the display list
155   - ///@param map is the mapping relationship between two networks
156   - ///@param colormap is the random generated color set for render
157   - void glRandColorCylinder1(GLuint &dlist1, std::vector<unsigned> map, std::vector<T> colormap, float sigma, float radius) {
  202 + ///render a T as a adjoint network of GT in transparancy(darkgreen, overlaid)
  203 + void glAdjointCylinder(float sigma, float radius) {
158 204  
159   - if (radius != sigma) // if render radius was changed by user, create a new display list
160   - glDeleteLists(dlist1, 1);
161   - if (!glIsList(dlist1)) { // if dlist1 isn't a display list, create it
162   - dlist1 = glGenLists(1); // generate a display list
163   - glNewList(dlist1, GL_COMPILE); // start a new display list
164   - for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network
165   - if (map[e] != unsigned(-1)) {
166   - glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]);
167   - for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge
168   - stim::circle<T> C1 = E[e].circ(p - 1);
169   - stim::circle<T> C2 = E[e].circ(p);
170   - C1.set_R(2*radius); // scale the circle to the same
171   - C2.set_R(2*radius);
172   - std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
173   - std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
174   - renderCylinder(Cp1, Cp2);
175   - }
176   - }
177   - else {
178   - glColor3f(1.0f, 1.0f, 1.0f); // white color for the un-mapping edges
179   - for (unsigned p = 1; p < E[e].size(); p++) { // for each point on that edge
180   - stim::circle<T> C1 = E[e].circ(p - 1);
181   - stim::circle<T> C2 = E[e].circ(p);
182   - C1.set_R(2*radius); // scale the circle to the same
183   - C2.set_R(2*radius);
184   - std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
185   - std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
186   - renderCylinder(Cp1, Cp2);
  205 + if (radius != sigma) // if render radius was changed by user, create a new display list
  206 + glDeleteLists(dlist+4, 1);
  207 + if (!glIsList(dlist+4)) {
  208 + glNewList(dlist+4, GL_COMPILE);
  209 + for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network
  210 + for (unsigned p = 1; p < E[e].size(); p++) { // for each point on that edge
  211 + stim::circle<T> C1 = E[e].circ(p - 1);
  212 + stim::circle<T> C2 = E[e].circ(p);
  213 + C1.set_R(2 * radius); // scale the circle to the same
  214 + C2.set_R(2 * radius);
  215 + std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
  216 + std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
  217 + glBegin(GL_QUAD_STRIP);
  218 + for (unsigned i = 0; i < Cp1.size(); i++) { // for every point on the circle(+1 means closing the circle)
  219 + glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]);
  220 + glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]);
187 221 }
188   - }
  222 + glEnd();
  223 + } // set the texture coordinate based on the specified magnitude index
189 224 }
190   - for (unsigned v = 0; v < V.size(); v++) {
191   - size_t num_edge = V[v].e[0].size() + V[v].e[1].size();
192   - if (num_edge > 1) { // if it is the joint vertex
193   - glColor3f(0.3, 0.3, 0.3); // gray color
194   - renderBall(V[v][0], V[v][1], V[v][2], 3*radius, 20);
  225 + for (unsigned n = 0; n < V.size(); n++) {
  226 + size_t num = V[n].e[0].size(); // store the number of outgoing edge of that vertex
  227 + if (num != 0) { // if it has outgoing edge
  228 + unsigned idx = V[n].e[0][0]; // find the index of first outgoing edge of that vertex
195 229 }
196   - else { // if it is the terminal vertex
197   - glColor3f(0.6, 0.6, 0.6); // more white gray
198   - renderBall(V[v][0], V[v][1], V[v][2], 3*radius, 20);
  230 + else {
  231 + unsigned idx = V[n].e[1][0]; // find the index of first incoming edge of that vertex
199 232 }
  233 + renderBall(V[n][0], V[n][1], V[n][2], 2 * radius, 20, 20);
200 234 }
201 235 glEndList();
202 236 }
203   - glCallList(dlist1);
  237 + glCallList(dlist+4);
204 238 }
205 239  
206   - ///render the T network cylinder as series of tubes
207   - ///@param dlist2 is the display list
  240 + ///render the network cylinder as series of tubes
  241 + ///@param I is a indicator: 0 -> GT, 1 -> T
208 242 ///@param map is the mapping relationship between two networks
209 243 ///@param colormap is the random generated color set for render
210   - void glRandColorCylinder2(GLuint &dlist2, std::vector<unsigned> map, std::vector<T> colormap, float sigma, float radius) {
  244 + void glRandColorCylinder(int I, std::vector<unsigned> map, std::vector<T> colormap, float sigma, float radius) {
211 245  
212   - if (radius != sigma) // if render radius was changed by user, create a new display list
213   - glDeleteLists(dlist2, 1);
214   - if (!glIsList(dlist2)) {
215   - dlist2 = glGenLists(1);
216   - glNewList(dlist2, GL_COMPILE);
217   - for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network
  246 + if (radius != sigma) // if render radius was changed by user, create a new display list
  247 + glDeleteLists(dlist+2, 1);
  248 + if (!glIsList(dlist+2)) { // if dlist isn't a display list, create it
  249 + glNewList(dlist+2, GL_COMPILE); // start a new display list
  250 + glEnable(GL_COLOR_MATERIAL);
  251 + glColorMaterial(GL_FRONT, GL_AMBIENT_AND_DIFFUSE);
  252 + for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network
218 253 if (map[e] != unsigned(-1)) {
219   - glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]);
220   - for (unsigned p = 1; p < E[e].size(); p++) { // for each point on that edge(except last one)
  254 + if(I == 0) // if it is to render GT
  255 + glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]);
  256 + else // if it is to render T
  257 + glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]);
  258 +
  259 + for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge
221 260 stim::circle<T> C1 = E[e].circ(p - 1);
222 261 stim::circle<T> C2 = E[e].circ(p);
223   - C1.set_R(2*radius); // scale the circle to the same
  262 + C1.set_R(2*radius); // scale the circle to the same
224 263 C2.set_R(2*radius);
225 264 std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
226 265 std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
... ... @@ -244,60 +283,35 @@ public:
244 283 size_t num_edge = V[v].e[0].size() + V[v].e[1].size();
245 284 if (num_edge > 1) { // if it is the joint vertex
246 285 glColor3f(0.3, 0.3, 0.3); // gray color
247   - renderBall(V[v][0], V[v][1], V[v][2], 3*radius, 20);
  286 + renderBall(V[v][0], V[v][1], V[v][2], 3*radius, 20, 20);
248 287 }
249 288 else { // if it is the terminal vertex
250 289 glColor3f(0.6, 0.6, 0.6); // more white gray
251   - renderBall(V[v][0], V[v][1], V[v][2], 3*radius, 20);
  290 + renderBall(V[v][0], V[v][1], V[v][2], 3*radius, 20, 20);
252 291 }
253 292 }
254 293 glEndList();
  294 + glDisable(GL_COLOR_MATERIAL);
255 295 }
256   - glCallList(dlist2);
  296 + glCallList(dlist+2);
257 297 }
258 298  
259   - /// Render the GT network centerline as a series of line strips in random different color
  299 + ///render the network centerline as a series of line strips in random different color
  300 + ///@param I is a indicator: 0 -> GT, 1 -> T
260 301 ///@param dlist1 is the display list
261 302 ///@param map is the mapping relationship between two networks
262 303 ///@param colormap is the random generated color set for render
263   - void glRandColorCenterline1(GLuint &dlist1, std::vector<unsigned> map, std::vector<T> colormap) {
  304 + void glRandColorCenterline(int I, GLuint &dlist1, std::vector<unsigned> map, std::vector<T> colormap) {
264 305 if (!glIsList(dlist1)) {
265 306 dlist1 = glGenLists(1);
266 307 glNewList(dlist1, GL_COMPILE);
267 308 for (unsigned e = 0; e < E.size(); e++) {
268   - if (map[e] != unsigned(-1)) { // if it has corresponding edge in another network
269   - glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]);
270   - glBegin(GL_LINE_STRIP);
271   - for (unsigned p = 0; p < E[e].size(); p++) {
272   - glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
273   - }
274   - glEnd();
275   - }
276   - else {
277   - glColor3f(1.0, 1.0, 1.0); // white color
278   - glBegin(GL_LINE_STRIP);
279   - for (unsigned p = 0; p < E[e].size(); p++) {
280   - glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
281   - }
282   - glEnd();
283   - }
284   - }
285   - glEndList();
286   - }
287   - glCallList(dlist1);
288   - }
  309 + if (map[e] != unsigned(-1)) { // if it has corresponding edge in another network
  310 + if (I == 0) // if it is to render GT
  311 + glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]);
  312 + else // if it is to render T
  313 + glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]);
289 314  
290   - /// Render the T network centerline as a series of line strips in random different color
291   - ///@param dlist2 is the display list
292   - ///@param map is the mapping relationship between two networks
293   - ///@param colormap is the random generated color set for render
294   - void glRandColorCenterline2(GLuint &dlist2, std::vector<unsigned> map, std::vector<T> colormap) {
295   - if (!glIsList(dlist2)) {
296   - dlist2 = glGenLists(1);
297   - glNewList(dlist2, GL_COMPILE);
298   - for (unsigned e = 0; e < E.size(); e++) {
299   - if (map[e] != unsigned(-1)) { // if it has corresponding edge in another network
300   - glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]);
301 315 glBegin(GL_LINE_STRIP);
302 316 for (unsigned p = 0; p < E[e].size(); p++) {
303 317 glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
... ... @@ -305,7 +319,7 @@ public:
305 319 glEnd();
306 320 }
307 321 else {
308   - glColor3f(1.0, 1.0, 1.0); // white color
  322 + glColor3f(1.0, 1.0, 1.0); // white color
309 323 glBegin(GL_LINE_STRIP);
310 324 for (unsigned p = 0; p < E[e].size(); p++) {
311 325 glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
... ... @@ -315,7 +329,7 @@ public:
315 329 }
316 330 glEndList();
317 331 }
318   - glCallList(dlist2);
  332 + glCallList(dlist1);
319 333 }
320 334  
321 335 //void glRandColorCenterlineGT(GLuint &dlist1, std::vector<unsigned> map, std::vector<T> colormap){
... ...
stim/visualization/swc.h
  Show/Hide comments   View file @18f207c
... ... @@ -86,6 +86,8 @@ namespace stim {
86 86 class swc {
87 87 public:
88 88 std::vector< typename swc_tree::swc_node<T> > node;
  89 + std::vector< std::vector<int> > E; // list of nodes
  90 +
89 91 swc() {}; // default constructor
90 92  
91 93 // load the swc as tree nodes
... ... @@ -170,11 +172,104 @@ namespace stim {
170 172 }
171 173 }
172 174  
  175 + // create a new edge
  176 + void create_up(std::vector<int>& edge,swc_tree::swc_node<T> cur_node, int target) {
  177 + while (cur_node.parent_idx != target) {
  178 + edge.push_back(cur_node.idx);
  179 + cur_node = node[cur_node.parent_idx - 1]; // move to next node
  180 + }
  181 + edge.push_back(cur_node.idx); // push back the start/end vertex of current edge
  182 +
  183 + std::reverse(edge.begin(), edge.end()); // follow the original flow direction
  184 + }
  185 + void create_down(std::vector<int>& edge, swc_tree::swc_node<T> cur_node, int j) {
  186 + while (cur_node.son_idx.size() != 0) {
  187 + edge.push_back(cur_node.idx);
  188 + if (cur_node.son_idx.size() > 1)
  189 + cur_node = node[cur_node.son_idx[j] - 1]; // move to next node
  190 + else
  191 + cur_node = node[cur_node.son_idx[0] - 1];
  192 + }
  193 + edge.push_back(cur_node.idx); // push back the start/end vertex of current edge
  194 + }
  195 +
  196 + // resample the tree-like SWC
  197 + void resample() {
  198 + unsigned n = node.size();
  199 +
  200 + std::vector<int> joint_node;
  201 + for (unsigned i = 1; i < n; i++) { // search all nodes(except the first one) to find joint nodes
  202 + if (node[i].son_idx.size() > 1) {
  203 + joint_node.push_back(node[i].idx); // store the original index
  204 + }
  205 + }
  206 +
  207 + std::vector<int> new_edge; // new edge in the network
  208 +
  209 + n = joint_node.size();
  210 +
  211 + for (unsigned i = 0; i < n; i++) { // for every joint nodes
  212 + std::vector<int> new_edge; // new edge in the network
  213 + swc_tree::swc_node<T> cur_node = node[joint_node[i] - 1]; // store current node
  214 +
  215 + // go up
  216 + swc_tree::swc_node<T> tmp_node = node[cur_node.parent_idx - 1];
  217 + while (tmp_node.parent_idx != -1 && tmp_node.son_idx.size() == 1) {
  218 + tmp_node = node[tmp_node.parent_idx - 1];
  219 + }
  220 + int target = tmp_node.parent_idx; // store the go-up target
  221 + create_up(new_edge, cur_node, target);
  222 + E.push_back(new_edge);
  223 + new_edge.clear();
  224 +
  225 + // go down
  226 + unsigned t = cur_node.son_idx.size();
  227 + for (unsigned j = 0; j < t; j++) { // for every son node
  228 + tmp_node = node[cur_node.son_idx[j] - 1]; // move down
  229 + while (tmp_node.son_idx.size() == 1) {
  230 + tmp_node = node[tmp_node.son_idx[0] - 1];
  231 + }
  232 + if (tmp_node.son_idx.size() == 0) {
  233 + create_down(new_edge, cur_node, j);
  234 + E.push_back(new_edge);
  235 + new_edge.clear();
  236 + }
  237 + }
  238 + }
  239 + }
  240 +
  241 + // get points in one edge
  242 + void get_points(int e, std::vector< stim::vec3<T> >& V) {
  243 + V.resize(E[e].size());
  244 + for (unsigned i = 0; i < E[e].size(); i++) {
  245 + unsigned id = E[e][i] - 1;
  246 +
  247 + V[i][0] = node[id].point[0];
  248 + V[i][1] = node[id].point[1];
  249 + V[i][2] = node[id].point[2];
  250 + }
  251 + }
  252 +
  253 + // get radius information in one edge
  254 + void get_radius(int e, std::vector<T>& radius) {
  255 + radius.resize(E[e].size());
  256 + for (unsigned i = 0; i < E[e].size(); i++) {
  257 + unsigned id = E[e][i] - 1;
  258 +
  259 + radius[i] = node[id].radius;
  260 + }
  261 + }
  262 +
173 263 // return the number of point in swc
174   - unsigned int numL() {
  264 + unsigned int numP() {
175 265 return node.size();
176 266 }
177 267  
  268 + // return the number of edges in swc after resample
  269 + unsigned int numE() {
  270 + return E.size();
  271 + }
  272 +
178 273  
179 274 };
180 275 } // end of namespace stim
... ...