1599ae65
Jiaming Guo
add swc.h and fun...
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
|
#ifndef STIM_SWC_H
#define STIM_SWC_H
#include <vector>
#include <fstream>
#include <sstream>
#include <iostream>
//STIM includes
#include <stim/math/vec3.h>
#include <stim/parser/parser.h>
namespace stim {
namespace swc_tree {
template <typename T>
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:
|
58674036
Jiaming Guo
fix minor error f...
|
23
|
int idx; // the index of current node start from 1(original index from the file)
|
1599ae65
Jiaming Guo
add swc.h and fun...
|
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
|
neuronal_type type; // the type of neuronal segmemt
stim::vec3<T> 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<int> 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<std::string> 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
|
58674036
Jiaming Guo
fix minor error f...
|
47
|
|
1599ae65
Jiaming Guo
add swc.h and fun...
|
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
|
// 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 <typename T>
class swc {
public:
|
ee17b90b
Jiaming Guo
make it compatibl...
|
88
|
std::vector< typename swc_tree::swc_node<T> > node;
|
2b28db3e
Jiaming Guo
better render
|
89
90
|
std::vector< std::vector<int> > E; // list of nodes
|
1599ae65
Jiaming Guo
add swc.h and fun...
|
91
92
93
94
95
96
97
98
99
100
|
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) {
|
da89bce9
Jiaming Guo
fixed splitting bugs
|
101
|
std::cerr << "STIM::SWC Error loading file" << filename << std::endl;
|
1599ae65
Jiaming Guo
add swc.h and fun...
|
102
103
104
105
106
107
|
exit(-1);
}
std::string line;
// skip comment
while (getline(infile, line)) {
|
58674036
Jiaming Guo
fix minor error f...
|
108
109
|
if ('#' == line[0] || line.empty()) // if it is comment line or empty line
continue; // keep read
|
1599ae65
Jiaming Guo
add swc.h and fun...
|
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
|
else
break;
}
unsigned int l = 0; // number of nodes
// get rid of the first/origin node
swc_tree::swc_node<T> 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<T> 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) {
|
da89bce9
Jiaming Guo
fixed splitting bugs
|
143
|
std::cerr << "STIM::SWC Error loading file" << filename << std::endl;
|
1599ae65
Jiaming Guo
add swc.h and fun...
|
144
145
146
147
148
|
exit(1);
}
std::string line;
while (getline(infile, line)) {
|
58674036
Jiaming Guo
fix minor error f...
|
149
|
if ('#' == line[0] || line.empty()) {
|
1599ae65
Jiaming Guo
add swc.h and fun...
|
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
|
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;
|
58674036
Jiaming Guo
fix minor error f...
|
170
|
int tmp_parent_idx = node[i].parent_idx - 1; // get the parent node loop idx of current node
|
1599ae65
Jiaming Guo
add swc.h and fun...
|
171
172
173
174
|
node[tmp_parent_idx].son_idx.push_back(i + 1); // son_idx stores the real idx = loop idx + 1
}
}
|
2b28db3e
Jiaming Guo
better render
|
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
|
// create a new edge
void create_up(std::vector<int>& edge,swc_tree::swc_node<T> 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<int>& edge, swc_tree::swc_node<T> 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() {
|
eceb6571
Jiaming Guo
fixed bug when sw...
|
198
|
unsigned nn = node.size();
|
2b28db3e
Jiaming Guo
better render
|
199
200
|
std::vector<int> joint_node;
|
eceb6571
Jiaming Guo
fixed bug when sw...
|
201
|
for (unsigned i = 1; i < nn; i++) { // search all nodes(except the first one) to find joint nodes
|
2b28db3e
Jiaming Guo
better render
|
202
203
204
205
206
207
208
|
if (node[i].son_idx.size() > 1) {
joint_node.push_back(node[i].idx); // store the original index
}
}
std::vector<int> new_edge; // new edge in the network
|
eceb6571
Jiaming Guo
fixed bug when sw...
|
209
|
unsigned n = joint_node.size();
|
2b28db3e
Jiaming Guo
better render
|
210
211
|
for (unsigned i = 0; i < n; i++) { // for every joint nodes
|
eceb6571
Jiaming Guo
fixed bug when sw...
|
212
|
|
2b28db3e
Jiaming Guo
better render
|
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
|
swc_tree::swc_node<T> cur_node = node[joint_node[i] - 1]; // store current node
// go up
swc_tree::swc_node<T> 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();
}
}
}
|
eceb6571
Jiaming Guo
fixed bug when sw...
|
239
240
241
242
243
244
245
|
if (n == 0) { // just one edge
new_edge.clear();
for (unsigned i = 0; i < nn; i++)
new_edge.push_back(node[i].idx);
E.push_back(new_edge);
}
|
2b28db3e
Jiaming Guo
better render
|
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
|
}
// get points in one edge
void get_points(int e, std::vector< stim::vec3<T> >& 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<T>& 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;
}
}
|
1599ae65
Jiaming Guo
add swc.h and fun...
|
270
|
// return the number of point in swc
|
2b28db3e
Jiaming Guo
better render
|
271
|
unsigned int numP() {
|
1599ae65
Jiaming Guo
add swc.h and fun...
|
272
273
274
|
return node.size();
}
|
2b28db3e
Jiaming Guo
better render
|
275
276
277
278
279
|
// return the number of edges in swc after resample
unsigned int numE() {
return E.size();
}
|
1599ae65
Jiaming Guo
add swc.h and fun...
|
280
281
282
283
284
|
};
} // end of namespace stim
#endif
|