121ccba7
pranathivemuri
network class and...
|
1
2
3
4
5
6
7
|
#ifndef STIM_FIBER_H
#define STIM_FIBER_H
#include <vector>
#include <ANN/ANN.h>
namespace stim{
|
3800c27f
David Mayerich
added code to sub...
|
8
|
|
121ccba7
pranathivemuri
network class and...
|
9
10
11
12
13
14
15
|
/** This class stores information about a single fiber represented as a set of geometric points
* between two branch or end points. This class is used as a fundamental component of the stim::network
* class to describe an interconnected (often biological) network.
*/
template<typename T>
class fiber{
|
7f27eafa
David Mayerich
simplified the st...
|
16
|
protected:
|
121ccba7
pranathivemuri
network class and...
|
17
18
19
20
|
unsigned int N; //number of points in the fiber
double **c; //centerline (array of double pointers)
T* r; // array of fiber radii
|
43eddb69
pranathivemuri
changed length's ...
|
21
|
ANNkd_tree* kdt; //kd-tree stores all points in the fiber for fast searching
|
121ccba7
pranathivemuri
network class and...
|
22
23
|
/// Initialize an empty fiber
|
964b4898
pranathivemuri
added function to...
|
24
25
|
void init()
{
|
121ccba7
pranathivemuri
network class and...
|
26
27
28
29
30
31
32
|
kdt = NULL;
c=NULL;
r=NULL;
N=0;
}
/// Initialize a fiber with N centerline points (all located at [0, 0, 0] with radius 0)
|
964b4898
pranathivemuri
added function to...
|
33
34
|
void init(unsigned int n)
{
|
121ccba7
pranathivemuri
network class and...
|
35
36
37
38
39
40
41
|
N = n; //set the number of points
kdt = NULL;
c = (double**) malloc(sizeof(double*) * N); //allocate the array pointer
for(unsigned int i = 0; i < N; i++) //allocate space for each point
c[i] = (double*) malloc(sizeof(double) * 3);
|
3800c27f
David Mayerich
added code to sub...
|
42
|
|
121ccba7
pranathivemuri
network class and...
|
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
|
r = (T*) malloc(sizeof(T) * N); //allocate space for the radii
}
/// Copies an existing fiber to the current fiber
/// @param cpy stores the new copy of the fiber
void copy( const stim::fiber<T>& cpy ){
///allocate space for the new fiber
init(cpy.N);
///copy the points
for(unsigned int i = 0; i < N; i++){
for(unsigned int d = 0; d < 3; d++) //for each dimension
c[i][d] = cpy.c[i][d]; //copy the coordinate
r[i] = cpy.r[i]; //copy the radius
}
gen_kdtree(); //generate the kd tree for the new fiber
}
|
3800c27f
David Mayerich
added code to sub...
|
65
|
/// generate a KD tree for points on fiber
|
964b4898
pranathivemuri
added function to...
|
66
67
68
|
void gen_kdtree()
{
int n_data = N; //create an array of data points
|
121ccba7
pranathivemuri
network class and...
|
69
|
ANNpointArray pts = (ANNpointArray)c; //cast the centerline list to an ANNpointArray
|
121ccba7
pranathivemuri
network class and...
|
70
71
72
|
kdt = new ANNkd_tree(pts, n_data, 3); //build a KD tree
}
|
964b4898
pranathivemuri
added function to...
|
73
|
/// find distance between two points
|
121ccba7
pranathivemuri
network class and...
|
74
75
|
double dist(double* p0, double* p1){
|
964b4898
pranathivemuri
added function to...
|
76
|
double sum = 0; // initialize variables
|
7f27eafa
David Mayerich
simplified the st...
|
77
|
float v;
|
964b4898
pranathivemuri
added function to...
|
78
79
|
for(unsigned int d = 0; d < 3; d++)
{
|
7f27eafa
David Mayerich
simplified the st...
|
80
81
82
83
|
v = p1[d] - p0[d];
sum +=v * v;
}
|
121ccba7
pranathivemuri
network class and...
|
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
|
return sqrt(sum);
}
/// This function retreives the index for the fiber point closest to q
/// @param q is a reference point used to find the closest point on the fiber center line
unsigned int ann( stim::vec<double> q ){
ANNidxArray idx = new ANNidx[1]; //variable used to hold the nearest point
ANNdistArray sq_dist = new ANNdist[1]; //variable used to hold the squared distance to the nearest point
kdt->annkSearch(q.data(), 1, idx, sq_dist); //search the KD tree for the nearest neighbor
return *idx;
}
|
7f27eafa
David Mayerich
simplified the st...
|
100
101
102
103
104
105
106
107
108
109
110
111
112
|
/// Returns a stim::vec representing the point at index i
/// @param i is an index of the desired centerline point
stim::vec<T> get_vec(unsigned i){
stim::vec<T> r;
r.resize(3);
r[0] = c[i][0];
r[1] = c[i][1];
r[2] = c[i][2];
return r;
}
|
121ccba7
pranathivemuri
network class and...
|
113
|
|
121ccba7
pranathivemuri
network class and...
|
114
115
116
117
118
119
120
121
122
123
124
125
|
public:
fiber(){
init();
}
/// Copy constructor
fiber(const stim::fiber<T> &obj){
copy(obj);
}
|
3800c27f
David Mayerich
added code to sub...
|
126
|
|
121ccba7
pranathivemuri
network class and...
|
127
128
129
130
131
132
|
//temp constructor for graph visualization
fiber(int n)
{
init(n);
}
|
7f27eafa
David Mayerich
simplified the st...
|
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
|
/// Constructor takes a list of stim::vec points, the radius at each point is set to zero
fiber(std::vector< stim::vec<T> > p){
init(p.size()); //initialize the fiber
//for each point, set the centerline position and radius
for(unsigned int i = 0; i < N; i++){
//set the centerline position
for(unsigned int d = 0; d < 3; d++)
c[i][d] = (double) p[i][d];
//set the radius
r[i] = 0;
}
//generate a kd tree
gen_kdtree();
}
|
121ccba7
pranathivemuri
network class and...
|
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
|
/// constructor takes a list of points and radii
fiber(std::vector< stim::vec< T > > pos, std::vector< T > radii){
init(pos.size()); //initialize the fiber
//for each point, set the centerline position and radius
for(unsigned int i = 0; i < N; i++){
//set the centerline position
for(unsigned int d = 0; d < 3; d++)
c[i][d] = (double) pos[i][d];
//set the radius
r[i] = radii[i];
}
//generate a kd tree
gen_kdtree();
}
|
121ccba7
pranathivemuri
network class and...
|
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
|
/// constructor takes an array of points and radii
// this function is used when the radii are represented as a stim::vec,
// since this may be easier when importing OBJs
fiber(std::vector< stim::vec<T> > pos, std::vector< stim::vec<T> > radii){
init(pos.size());
//for each point, set the position and radius
for(unsigned int i = 0; i < N; i++){
//at(i) = (double*)malloc(sizeof(double) * 3);
for(unsigned int d = 0; d < 3; d++)
c[i][d] = (double) pos[i][d];
r[i] = radii[i][(unsigned int)0];
}
gen_kdtree();
}
/// Assignment operation
fiber& operator=(const fiber &rhs){
if(this == &rhs) return *this; //test for and handle self-assignment
copy(rhs);
|
121ccba7
pranathivemuri
network class and...
|
196
197
198
|
}
/// Calculate the length of the fiber and return it.
|
43eddb69
pranathivemuri
changed length's ...
|
199
|
double length(){
|
121ccba7
pranathivemuri
network class and...
|
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
|
double* p0;
double *p1;
double l = 0; //initialize the length to zero
//for each point
//typename std::list< point<T> >::iterator i; //create a point iterator
for(unsigned int i = 0; i < N; i++){ //for each point in the fiber
if(i == 0) //if this is the first point, just store it
p1 = c[0];
else{ //if this is any other point
p0 = p1; //shift p1->p0
p1 = c[i]; //set p1 to the new point
l += dist(p0, p1); //add the length of p1 - p0 to the running sum
}
}
|
43eddb69
pranathivemuri
changed length's ...
|
218
|
return l; //return the length
|
121ccba7
pranathivemuri
network class and...
|
219
|
}
|
3800c27f
David Mayerich
added code to sub...
|
220
|
|
121ccba7
pranathivemuri
network class and...
|
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
|
/// Calculates the length and average radius of the fiber
/// @param length is filled with the fiber length
T radius(T& length){
double* p0; //temporary variables to store point positions
double* p1;
T r0, r1; //temporary variables to store radii at points
double l;
T r_mean; //temporary variable to store the length and average radius of a fiber segment
double 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<T> >::iterator i; //create a point iterator
for(unsigned int i = 0; i < N; i++){ //for each point in the fiber
if(i == 0){ //if this is the first point, just store it
p1 = c[0];
r1 = r[0];
}
else{ //if this is any other point
p0 = p1; //shift p1->p0 and r1->r0
r0 = r1;
p1 = c[i]; //set p1 to the new point
r1 = r[i];
l = dist(p0, p1); //calculate the length of the p0-p1 segment
r_mean = (r0 + r1) / 2; //calculate the average radius of the segment
radius_sum += r_mean * (T) 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 (T)(radius_sum / length); //return the average radius of the fiber
}
T average_radius()
{
T r_sum = 0.;
for(unsigned int i = 0; i < N; i++)
{
|
3800c27f
David Mayerich
added code to sub...
|
269
|
r_sum = r_sum + r[i];
|
121ccba7
pranathivemuri
network class and...
|
270
271
272
273
274
275
276
277
278
279
280
281
282
283
|
}
return r_sum/((T) N);
}
/// Calculates the average radius of the fiber
T radius(){
T length;
return radius(length);
}
/// Returns the radius at index idx.
T radius(int idx){
return r[idx];
}
|
121ccba7
pranathivemuri
network class and...
|
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
|
/// Return the point on the fiber closest to q
/// @param q is the query point used to locate the nearest point on the fiber centerline
stim::vec<T> nearest(stim::vec<T> q){
stim::vec<double> temp( (double) q[0], (double) q[1], (double) q[2]);
unsigned int idx = ann(temp); //determine the index of the nearest neighbor
return stim::vec<T>((T) c[idx][0], (T) c[idx][1], (T) c[idx][2]); //return the nearest centerline point
}
/// Return the point index on the fiber closest to q
/// @param q is the query point used to locate the nearest point on the fiber centerline
unsigned int nearest_idx(stim::vec<T> q){
|
4df9ec0e
pranathivemuri
compare networks ...
|
300
|
stim::vec<double> temp((double) q[0], (double) q[1], (double) q[2]);
|
121ccba7
pranathivemuri
network class and...
|
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
|
unsigned int idx = ann(temp); //determine the index of the nearest neighbor
return idx; //return the nearest centerline point index
}
/// Returns the fiber centerline as an array of stim::vec points
std::vector< stim::vec<T> > centerline(){
//create an array of stim vectors
std::vector< stim::vec<T> > pts(N);
//cast each point to a stim::vec, keeping only the position information
for(unsigned int i = 0; i < N; i++)
pts[i] = stim::vec<T>((T) c[i][0], (T) c[i][1], (T) c[i][2]);
//return the centerline array
return pts;
}
/// Returns the fiber centerline magnitudes as an array of stim::vec points
std::vector< stim::vec<T> > centerlinemag(){
//create an array of stim vectors
std::vector< stim::vec<T> > pts(N);
//cast each point to a stim::vec, keeping only the position information
for(unsigned int i = 0; i < N; i++)
pts[i] = stim::vec<T>(r[i], r[i]);;
//return the centerline array
return pts;
}
/// Split the fiber at the specified index. If the index is an end point, only one fiber is returned
std::vector< stim::fiber<T> > split(unsigned int idx){
std::vector< stim::fiber<T> > fl; //create an array to store up to two fibers
//if the index is an end point, only the existing fiber is returned
if(idx == 0 || idx == N-1){
fl.resize(1); //set the size of the fiber to 1
fl[0] = *this; //copy the current fiber
}
//if the index is not an end point
else{
unsigned int N1 = idx + 1; //calculate the size of both fibers
unsigned int N2 = N - idx;
fl.resize(2); //set the array size to 2
fl[0].init(N1); //set the size of each fiber
fl[1].init(N2);
//copy both halves of the fiber
unsigned int i, d;
//first half
for(i = 0; i < N1; i++){ //for each centerline point
for(d = 0; d < 3; d++)
fl[0].c[i][d] = c[i][d]; //copy each coordinate
fl[0].r[i] = r[i]; //copy the corresponding radius
}
//second half
for(i = 0; i < N2; i++){
for(d = 0; d < 3; d++)
fl[1].c[i][d] = c[idx + i][d];
fl[1].r[i] = r[idx + i];
}
}
return fl; //return the array
}
/// Calculates the set of fibers resulting from a connection between the current fiber and a fiber f
/// @param f is the fiber that will be connected to the current fiber
std::vector< stim::fiber<T> > connect( stim::fiber<T> &f, double dist){
double min_dist;
unsigned int idx0, idx1;
//go through each point in the query fiber, looking for the indices for the closest points
for(unsigned int i = 0; i < f.n_pts(); i++){
//Run through all points and find the index with the closest point, then partition the fiber and return two fibers.
}
}
/// Outputs the fiber as a string
std::string str(){
std::stringstream ss;
//create an iterator for the point list
//typename std::list< point<T> >::iterator i;
for(unsigned int i = 0; i < N; i++){
ss<<" [ ";
for(unsigned int d = 0; d < 3; d++){
ss<<c[i][d]<<" ";
}
ss<<"] r = "<<r[i]<<std::endl;
}
return ss.str();
}
|
121ccba7
pranathivemuri
network class and...
|
413
|
/// Returns the number of centerline points in the fiber
|
797f8ab0
David Mayerich
cleaned up the st...
|
414
|
unsigned int size(){
|
121ccba7
pranathivemuri
network class and...
|
415
416
|
return N;
}
|
7f27eafa
David Mayerich
simplified the st...
|
417
|
|
4df9ec0e
pranathivemuri
compare networks ...
|
418
|
|
7f27eafa
David Mayerich
simplified the st...
|
419
420
421
422
|
/// Bracket operator returns the element at index i
/// @param i is the index of the element to be returned as a stim::vec
stim::vec<T> operator[](unsigned i){
|
3800c27f
David Mayerich
added code to sub...
|
423
|
return get_vec(i);
|
7f27eafa
David Mayerich
simplified the st...
|
424
425
426
427
428
429
|
}
/// Back method returns the last point in the fiber
stim::vec<T> back(){
return get_vec(N-1);
}
|
964b4898
pranathivemuri
added function to...
|
430
|
////resample a fiber in the network
|
797f8ab0
David Mayerich
cleaned up the st...
|
431
|
stim::fiber<T> resample(T spacing)
|
964b4898
pranathivemuri
added function to...
|
432
|
{
|
5bccf89d
David Mayerich
fixed some errors...
|
433
|
std::cout<<"fiber::resample()"<<std::endl;
|
7c43b7f9
pranathivemuri
return resampled ...
|
434
|
|
964b4898
pranathivemuri
added function to...
|
435
436
|
std::vector<T> v(3); //v-direction vector of the segment
stim::vec<T> p(3); //- intermediate point to be added
|
4df9ec0e
pranathivemuri
compare networks ...
|
437
438
|
stim::vec<T> p1(3); // p1 - starting point of an segment on the fiber,
stim::vec<T> p2(3); // p2 - ending point,
|
964b4898
pranathivemuri
added function to...
|
439
|
double sum=0; //distance summation
|
7c43b7f9
pranathivemuri
return resampled ...
|
440
|
std::vector<stim::vec<T> > fiberPositions = centerline();
|
964b4898
pranathivemuri
added function to...
|
441
442
|
std::vector<stim::vec<T> > newPointList; // initialize list of new resampled points on the fiber
// for each point on the centerline (skip if it is the last point on centerline)
|
797f8ab0
David Mayerich
cleaned up the st...
|
443
|
//unsigned int N = fiberPositions.size(); // number of points on the fiber
|
964b4898
pranathivemuri
added function to...
|
444
445
|
for(unsigned int f=0; f< N-1; f++)
{
|
4df9ec0e
pranathivemuri
compare networks ...
|
446
447
448
449
|
p1 = fiberPositions[f]; p2 = fiberPositions[f + 1]; v = p2 - p1;
for(unsigned int d = 0; d < 3; d++){
sum +=v[d] * v[d];} //length of segment-distance between starting and ending point
|
964b4898
pranathivemuri
added function to...
|
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
|
T lengthSegment = sqrt(sum); //find Length of the segment as distance between the starting and ending points of the segment
if(lengthSegment >= spacing) // if length of the segment is greater than standard deviation resample
{
// repeat resampling until accumulated stepsize is equsl to length of the segment
for(T step=0.0; step<lengthSegment; step+=spacing)
{
// calculate the resampled point by travelling step size in the direction of normalized gradient vector
for(unsigned int i=0; i<3;i++)
{
p[i] = p1[i] + v[i]*(step/lengthSegment);
} //for each dimension
// add this resampled points to the new fiber list
newPointList.push_back(p);
}
}
else // length of the segment is now less than standard deviation, push the ending point of the segment and proceed to the next point in the fiber
newPointList.push_back(fiberPositions[f+1]);
}
newPointList.push_back(fiberPositions[N-1]); //add the last point on the fiber to the new fiber list
|
7c43b7f9
pranathivemuri
return resampled ...
|
471
472
|
fiber newFiber(newPointList);
return newFiber;
|
964b4898
pranathivemuri
added function to...
|
473
|
}
|
7f27eafa
David Mayerich
simplified the st...
|
474
|
|
121ccba7
pranathivemuri
network class and...
|
475
476
477
478
479
480
481
482
483
|
};
} //end namespace stim
#endif
|