df036f4c
 David Mayerich
added the network...
|
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
|
#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 > 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 = P.begin(); i != P.end(); i++){ //for each point in the fiber
if(i == P.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 = P.begin(); i != P.end(); i++){ //for each point in the fiber
if(i == P.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
return radius_sum / length; //return the average radius of the fiber
}
std::string str(){
std::stringstream ss;
//create an iterator for the point list
typename std::list<point>::iterator i;
for(i = P.begin(); i != P.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();
}
};
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
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->P.push_back(p); //push the point onto the current fiber
}
}
} //end load()
/// 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
|
Added #include <s...