network.h
12.6 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
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
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
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
143
144
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
203
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
262
263
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
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
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
414
415
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
#ifndef STIM_NETWORK_H
#define STIM_NETWORK_H
#include <stim/math/vector.h>
#include <stim/visualization/obj.h>
#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{
//-------------------HELPER CLASSES-----------------------
/// 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>
class fiber : public std::list< point >{
using std::list< point >::begin;
using std::list< point >::end;
using std::list< point >::size;
public:
//std::list< point > P; //geometric point positions
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
for(i = begin(); i != end(); i++){ //for each point in the fiber
if(i == begin()) //if this is the first point, just store it
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
for(i = begin(); i != end(); i++){ //for each point in the fiber
if(i == begin()){ //if this is the first point, just store it
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
//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
}
std::vector< stim::vec<T> > geometry(){
std::vector< stim::vec<T> > result; //create an array to store the fiber geometry
result.resize( size() ); //pre-allocate the array
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
for(p = begin(); p != end(); p++){
result[pi] = *p;
pi++;
}
return result; //return the geometry array
}
std::string str(){
std::stringstream ss;
//create an iterator for the point list
typename std::list<point>::iterator i;
for(i = begin(); i != end(); i++){
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();
}
};
//---------------NETWORK CLASS-----------------------------
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
/// @param object is the object file to be used as the basis for the network
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
f->push_back(p); //push the point onto the current fiber
}
}
} //end load()
/// 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)
result[vi] = ni->out.front()->front();
else if(ni->in.size() != 0)
result[vi] = ni->in.front()->back();
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;
}
/// 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