df036f4c
David Mayerich
added the network...
|
1
2
3
|
#ifndef STIM_NETWORK_H
#define STIM_NETWORK_H
|
c61fd046
David Mayerich
fixed vessel mode...
|
4
|
#include <stim/math/vector.h>
|
5687cf1b
Cherub P. Harder
Added #include <s...
|
5
|
#include <stim/visualization/obj.h>
|
df036f4c
David Mayerich
added the network...
|
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
|
#include <list>
#include <ANN/ANN.h>
namespace stim{
/** This class provides an interface for dealing with biological networks.
* It takes the following aspects into account:
* 1) Network geometry and centerlines
* 2) Network connectivity (a graph structure can be extracted)
* 3) Network surface structure (the surface is represented as a triangular mesh and referenced to the centerline)
*/
template<typename T>
class network{
|
e376212f
David Mayerich
added support for...
|
21
|
//-------------------HELPER CLASSES-----------------------
|
df036f4c
David Mayerich
added the network...
|
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
|
/// Stores information about a geometric point on the network centerline (including point position and radius)
//template<typename T>
class point : public stim::vec<T>{
public:
T r;
point() : stim::vec<T>(){}
//casting constructor
point(stim::vec<T> rhs) : stim::vec<T>(rhs){}
};
//template<typename T>
class t_node;
class fiber;
//create typedefs for the iterators to simplify the network code
typedef typename std::list< fiber >::iterator fiber_i;
typedef typename std::list< t_node >::iterator t_node_i;
/// Stores information about a single capillary (a length of vessel between two branch or end points)
//template<typename T>
|
cbce396e
David Mayerich
added support for...
|
45
46
47
48
49
50
|
class fiber : public std::list< point >{
using std::list< point >::begin;
using std::list< point >::end;
using std::list< point >::size;
|
df036f4c
David Mayerich
added the network...
|
51
52
|
public:
|
cbce396e
David Mayerich
added support for...
|
53
|
//std::list< point > P; //geometric point positions
|
df036f4c
David Mayerich
added the network...
|
54
55
56
57
58
59
60
61
62
63
64
65
66
67
|
typename std::list< t_node >::iterator n[2]; //indices to terminal nodes
unsigned int id;
public:
/// Calculate the length of the fiber and return it.
T length(){
point p0, p1;
T l = 0; //initialize the length to zero
//for each point
typename std::list< point >::iterator i; //create a point iterator
|
cbce396e
David Mayerich
added support for...
|
68
|
for(i = begin(); i != end(); i++){ //for each point in the fiber
|
df036f4c
David Mayerich
added the network...
|
69
|
|
cbce396e
David Mayerich
added support for...
|
70
|
if(i == begin()) //if this is the first point, just store it
|
df036f4c
David Mayerich
added the network...
|
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
|
p1 = *i;
else{ //if this is any other point
p0 = p1; //shift p1->p0
p1 = *i; //set p1 to the new point
l += (p1 - p0).len(); //add the length of p1 - p0 to the running sum
}
}
return l; //return the length
}
T radius(T& length){
point p0, p1; //temporary variables to store point positions
T r0, r1; //temporary variables to store radii at points
T l, r; //temporary variable to store the length and average radius of a fiber segment
T length_sum = 0; //initialize the length to zero
T radius_sum = 0; //initialize the radius sum to zero
//for each point
typename std::list< point >::iterator i; //create a point iterator
|
cbce396e
David Mayerich
added support for...
|
92
|
for(i = begin(); i != end(); i++){ //for each point in the fiber
|
df036f4c
David Mayerich
added the network...
|
93
|
|
cbce396e
David Mayerich
added support for...
|
94
|
if(i == begin()){ //if this is the first point, just store it
|
df036f4c
David Mayerich
added the network...
|
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
|
p1 = *i;
r1 = i->r;
}
else{ //if this is any other point
p0 = p1; //shift p1->p0 and r1->r0
r0 = r1;
p1 = *i; //set p1 to the new point
r1 = i->r; //and r1
l = (p1 - p0).len(); //calculate the length of the p0-p1 segment
r = (r0 + r1) / 2; //calculate the average radius of the segment
radius_sum += r * l; //add the radius scaled by the length to a running sum
length_sum += l; //add the length of p1 - p0 to the running sum
}
}
length = length_sum; //store the total length
|
fa737592
David Mayerich
fixed a bug that ...
|
113
114
115
116
117
118
|
//if the total length is zero, store a radius of zero
if(length == 0)
return 0;
else
return radius_sum / length; //return the average radius of the fiber
|
df036f4c
David Mayerich
added the network...
|
119
120
|
}
|
e376212f
David Mayerich
added support for...
|
121
122
123
|
std::vector< stim::vec<T> > geometry(){
std::vector< stim::vec<T> > result; //create an array to store the fiber geometry
|
cbce396e
David Mayerich
added support for...
|
124
|
result.resize( size() ); //pre-allocate the array
|
e376212f
David Mayerich
added support for...
|
125
126
127
128
129
|
typename std::list< point >::iterator p; //create a list iterator
unsigned int pi = 0; //create an index into the result array
//for each geometric point on the fiber centerline
|
cbce396e
David Mayerich
added support for...
|
130
|
for(p = begin(); p != end(); p++){
|
e376212f
David Mayerich
added support for...
|
131
132
133
134
135
136
137
138
|
result[pi] = *p;
pi++;
}
return result; //return the geometry array
}
|
df036f4c
David Mayerich
added the network...
|
139
140
141
142
143
|
std::string str(){
std::stringstream ss;
//create an iterator for the point list
typename std::list<point>::iterator i;
|
cbce396e
David Mayerich
added support for...
|
144
|
for(i = begin(); i != end(); i++){
|
df036f4c
David Mayerich
added the network...
|
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
|
ss<<i->str()<<" r = "<<i->r<<std::endl;
}
return ss.str();
}
};
/// Terminal node for a capillary. This is analogous to a graph vertex and contains a list of edge indices.
//template<typename T>
class t_node{
public:
unsigned int id;
//lists of edge indices for capillaries
//the "in" and "out" just indicate how the geometry is defined:
// edges in the "in" list are geometrically oriented such that the terminal node is last
// edges in the "out" list are geometrically oriented such that the terminal node is first
std::list< fiber_i > in; //edge indices for incoming capillaries
std::list< fiber_i > out; //edge indices for outgoing capillaries
std::string str(){
std::stringstream ss;
ss<<id<<": "; //output the node ID
//output the IDs for both lists
typename std::list< fiber_i >::iterator f;
for(f = in.begin(); f != in.end(); f++){
if(f != in.begin())
ss<<", ";
ss<<(*f)->n[0]->id;
}
//if there are nodes in both lists, separate them by a comma
if(out.size() > 0 && in.size() > 0)
ss<<", ";
for(f = out.begin(); f != out.end(); f++){
if(f != out.begin())
ss<<", ";
ss<<(*f)->n[1]->id;
}
return ss.str();
}
};
|
e376212f
David Mayerich
added support for...
|
203
|
//---------------NETWORK CLASS-----------------------------
|
df036f4c
David Mayerich
added the network...
|
204
205
206
207
208
209
210
211
212
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
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
|
protected:
//list of terminal nodes
std::list<t_node> N;
//list of fibers
std::list<fiber> F;
/// Sets a unique ID for each terminal node and fiber
void set_names(){
unsigned int i;
i = 0;
for(t_node_i ti = N.begin(); ti != N.end(); ti++)
ti->id = i++;
i = 0;
for(fiber_i fi = F.begin(); fi != F.end(); fi++)
fi->id = i++;
}
public:
std::string str(){
//assign names to elements of the network
set_names();
//create a stringstream for output
std::stringstream ss;
//output the nodes
ss<<"Nodes----------------------------"<<std::endl;
for(t_node_i i = N.begin(); i != N.end(); i++){
ss<<i->str()<<std::endl;
}
//output the fibers
ss<<std::endl<<"Fibers---------------------------"<<std::endl;
T length, radius;
//output every fiber
for(fiber_i f = F.begin(); f != F.end(); f++){
//calculate the length and average radius
radius = f->radius(length);
//output the IDs of the terminal nodes
ss<<f->n[0]->id<<" -- "<<f->n[1]->id<<": length = "<<length<<", average radius = "<<radius<<std::endl;
}
return ss.str();
}
/// Load a network from an OBJ object
|
e376212f
David Mayerich
added support for...
|
262
263
|
/// @param object is the object file to be used as the basis for the network
|
df036f4c
David Mayerich
added the network...
|
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
|
void load( stim::obj<T> object){
//get the number of vertices in the object
unsigned int nV = object.numV();
//allocate an array of pointers to nodes, which will be used to preserve connectivity
//initiate all values to T.end()
std::vector< t_node_i > node_hash(nV, N.end());
unsigned int nL = object.numL(); //get the number of lines in the OBJ
//for each line in the OBJ structure
for(unsigned int li = 0; li < nL; li++){
F.push_back(fiber()); //push a new fiber onto the fiber list
fiber_i f = --(F.end()); //get an iterator to the new fiber
//----------Handle the terminating nodes for the fiber
//get the indices of the line vertices
std::vector< unsigned int > Li = object.getL_Vi(li);
unsigned int i0 = Li.front() - 1;
unsigned int i1 = Li.back() - 1;
//deal with the first end point of the capillary
if(node_hash[i0] != N.end()){ //if the node has been used before
(*f).n[0] = node_hash[i0]; //assign the node to the new capillary
(*node_hash[i0]).out.push_back(f); //add an out pointer to the existing node
}
else{ //otherwise
N.push_back(t_node()); //create a new node and add it to the node list
t_node_i t = --(N.end()); //get an iterator to the new node
node_hash[i0] = t; //add a pointer to the new node to the hash list
(*f).n[0] = t; //add a pointer to the new node to the capillary
(*t).out.push_back(f); //add a pointer to the capillary to the new node
}
//deal with the last end point of the capillary
if(node_hash[i1] != N.end()){
(*f).n[1] = node_hash[i1];
(*node_hash[i1]).in.push_back(f);
}
else{
N.push_back(t_node());
t_node_i t = --(N.end());
node_hash[i1] = t; //add the new node to the hash list
(*f).n[1] = t;
(*t).in.push_back(f);
}
//-------------Handle the geometric points for the fiber
std::vector< vec<T> > L = object.getL_V(li);
std::vector< vec<T> > R = object.getL_VT(li);
unsigned int nP = L.size(); //get the number of geometric points in the fiber
//for each vertex in the fiber
for(unsigned int pi = 0; pi < nP; pi++){
point p = (point)L[pi]; //move the geometric coordinates into a point structure
p.r = R[pi][0]; //store the radius
|
cbce396e
David Mayerich
added support for...
|
324
|
f->push_back(p); //push the point onto the current fiber
|
df036f4c
David Mayerich
added the network...
|
325
326
327
328
329
|
}
}
} //end load()
|
e376212f
David Mayerich
added support for...
|
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
|
/// Returns an array of node positions
std::vector< stim::vec<T> > get_node_positions(){
std::vector< stim::vec<T> > result; //create an array to store the result
result.resize(N.size()); //set the array size
t_node_i ni; //create a terminal node iterator
unsigned int vi = 0; //vertex index into the result array
//for every terminal node
for(ni = N.begin(); ni != N.end(); ni++){
//create a vector based on the node position
//if the number of outgoing nodes is nonzero
if(ni->out.size() != 0)
|
cbce396e
David Mayerich
added support for...
|
346
|
result[vi] = ni->out.front()->front();
|
e376212f
David Mayerich
added support for...
|
347
|
else if(ni->in.size() != 0)
|
cbce396e
David Mayerich
added support for...
|
348
|
result[vi] = ni->in.front()->back();
|
e376212f
David Mayerich
added support for...
|
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
|
vi++; //increment the array index
}
//return the resulting array
return result;
}
std::vector< stim::vec<T> > get_fiber_geometry( fiber_i f ){
return f->geometry();
}
/// Generate an OBJ file from the network
stim::obj<T> obj(){
//create an OBJ object
stim::obj<T> object;
//name the nodes
set_names();
//retrieve a list of terminal node positions
std::vector< stim::vec<T> > node_pos = get_node_positions();
//add the nodes to the obj file
object.addV(node_pos);
//counter for vertex indices in the object class
unsigned int nP;
//for each fiber
fiber_i fi; //create a fiber iterator
for(fi = F.begin(); fi != F.end(); fi++){
//get an array of fiber points
std::vector< stim::vec<T> > fiber_p = get_fiber_geometry(fi);
//create a subset of this array
typename std::vector< stim::vec<T> >::iterator start = fiber_p.begin() + 1;
typename std::vector< stim::vec<T> >::iterator end = fiber_p.end() - 1;
typename std::vector< stim::vec<T> > fiber_subset(start, end);
//add this subset to the geometry object
nP = object.addV(fiber_subset);
//create an array to hold vertex indices for a line
std::vector<unsigned int> line;
line.resize(fiber_p.size());
//add the terminal nodes to the line list (make sure to add 1 to make them compatible with the OBJ)
line[0] = fi->n[0]->id + 1;
line[line.size() - 1] = fi->n[1]->id + 1;
//add the intermediate vertex indices to the line array
for(unsigned int i = 0; i < fiber_subset.size(); i++){
line[1 + i] = nP + i;
}
//add the line list to the object class
object.addLine(line);
}
return object;
|
e376212f
David Mayerich
added support for...
|
414
415
|
}
|
df036f4c
David Mayerich
added the network...
|
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
|
/// This function returns the information necessary for a simple graph-based physical (ex. fluid) simulation.
/// @param n0 is a array which will contain the list of source nodes
/// @param n1 is a array which will contain the list of destination nodes
/// @param length is a array containing the lengths of fibers in the network
/// @param radius is a array containing the average radii of fibers in the network
void build_simgraph(std::vector<unsigned int>& n0, std::vector<unsigned int>& n1, std::vector<T>& length, std::vector<T>& radius){
//determine the number of fibers in the network
unsigned int nF = F.size();
//allocate the necessary space to store the fiber information
n0.resize(nF);
n1.resize(nF);
length.resize(nF);
radius.resize(nF);
//assign names (identifiers) to the network components
set_names();
//fill the arrays
unsigned int i = 0;
T l, r;
for(fiber_i f = F.begin(); f != F.end(); f++){
n0[i] = f->n[0]->id; //get the identifiers for the first and second nodes for the current fiber
n1[i] = f->n[1]->id;
r = f->radius(l); //get the length and radius of the capillary (calculated at the same time)
radius[i] = r; //store the radius in the output array
length[i] = l; //store the length in the output array
i++; //increment the array index
}
}
};
}; //end namespace stim
#endif
|