a9ed210f
Pavel Govyadinov
added centerline
|
1
2
3
4
5
|
#ifndef STIM_CENTERLINE_H
#define STIM_CENTERLINE_H
#include <vector>
#include <stim/math/vec3.h>
|
c01958ad
Pavel Govyadinov
merged the two st...
|
6
|
//#include <ANN/ANN.h>
|
a9ed210f
Pavel Govyadinov
added centerline
|
7
8
9
10
11
12
13
14
|
namespace stim{
/** 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>
|
c01958ad
Pavel Govyadinov
merged the two st...
|
15
|
class centerline{
|
a9ed210f
Pavel Govyadinov
added centerline
|
16
17
|
protected:
|
c01958ad
Pavel Govyadinov
merged the two st...
|
18
19
20
21
22
23
24
25
26
27
28
|
unsigned int N; //number of points in the fiber
double **c; //centerline (array of double pointers)
// ANNkd_tree* kdt; //kd-tree stores all points in the fiber for fast searching
/// Initialize an empty fiber
void init()
{
N=0;
c=NULL;
// kdt = NULL;
}
|
a9ed210f
Pavel Govyadinov
added centerline
|
29
|
|
c01958ad
Pavel Govyadinov
merged the two st...
|
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
|
/// Initialize a fiber with N centerline points (all located at [0, 0, 0] with radius 0)
void init(unsigned int n)
{
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);
}
/// Copies an existing fiber to the current fiber
/// @param cpy stores the new copy of the fiber
void copy( const stim::centerline<T>& cpy, bool kd = 0){
///allocate space for the new fiber
init(cpy.N);
|
a9ed210f
Pavel Govyadinov
added centerline
|
49
|
|
c01958ad
Pavel Govyadinov
merged the two st...
|
50
51
52
53
|
///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
|
a9ed210f
Pavel Govyadinov
added centerline
|
54
|
}
|
c01958ad
Pavel Govyadinov
merged the two st...
|
55
56
|
// if(kd)
// gen_kdtree(); //generate the kd tree for the new fiber
|
a9ed210f
Pavel Govyadinov
added centerline
|
57
58
|
}
|
c01958ad
Pavel Govyadinov
merged the two st...
|
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
|
/// generate a KD tree for points on fiber
// void gen_kdtree()
// {
// int n_data = N; //create an array of data points
// ANNpointArray pts = (ANNpointArray)c; //cast the centerline list to an ANNpointArray
// kdt = new ANNkd_tree(pts, n_data, 3); //build a KD tree
// }
/// find distance between two points
double dist(double* p0, double* p1){
double sum = 0; // initialize variables
float v;
for(unsigned int d = 0; d < 3; d++)
{
v = p1[d] - p0[d];
sum +=v * v;
}
return sqrt(sum);
|
c72184d1
Jiaming Guo
add many function...
|
78
|
}
|
a9ed210f
Pavel Govyadinov
added centerline
|
79
|
|
c01958ad
Pavel Govyadinov
merged the two st...
|
80
81
82
83
84
85
86
87
88
89
|
/// 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;
// }
|
a9ed210f
Pavel Govyadinov
added centerline
|
90
91
92
93
|
/// 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){
|
c01958ad
Pavel Govyadinov
merged the two st...
|
94
95
96
97
98
|
stim::vec3<T> r;
r.resize(3);
r[0] = c[i][0];
r[1] = c[i][1];
r[2] = c[i][2];
|
a9ed210f
Pavel Govyadinov
added centerline
|
99
|
|
c01958ad
Pavel Govyadinov
merged the two st...
|
100
|
return r;
|
a9ed210f
Pavel Govyadinov
added centerline
|
101
102
103
104
105
|
}
public:
|
c01958ad
Pavel Govyadinov
merged the two st...
|
106
|
centerline(){
|
c72184d1
Jiaming Guo
add many function...
|
107
|
init();
|
a9ed210f
Pavel Govyadinov
added centerline
|
108
|
}
|
c01958ad
Pavel Govyadinov
merged the two st...
|
109
110
111
112
|
/// Copy constructor
centerline(const stim::centerline<T> &obj){
copy(obj);
|
a9ed210f
Pavel Govyadinov
added centerline
|
113
114
|
}
|
c01958ad
Pavel Govyadinov
merged the two st...
|
115
116
117
118
|
//temp constructor for graph visualization
centerline(int n)
{
init(n);
|
a9ed210f
Pavel Govyadinov
added centerline
|
119
120
|
}
|
c01958ad
Pavel Govyadinov
merged the two st...
|
121
122
123
124
125
126
|
/// Constructor takes a list of stim::vec points, the radius at each point is set to zero
centerline(std::vector< stim::vec<T> > p, bool kd = 0){
init(p.size()); //initialize the fiber
//for each point, set the centerline position and radius
for(unsigned int i = 0; i < N; i++){
|
a9ed210f
Pavel Govyadinov
added centerline
|
127
|
|
c01958ad
Pavel Govyadinov
merged the two st...
|
128
129
130
131
132
133
134
135
136
|
//set the centerline position
for(unsigned int d = 0; d < 3; d++)
c[i][d] = (double) p[i][d];
//set the radius
}
//generate a kd tree
// if(kd)
// gen_kdtree();
|
a9ed210f
Pavel Govyadinov
added centerline
|
137
138
|
}
|
c01958ad
Pavel Govyadinov
merged the two st...
|
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
|
/// constructor takes a list of points
centerline(std::vector< stim::vec3< T > > pos, bool kd = 0){
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
}
//generate a kd tree
//if(kd)
// gen_kdtree();
}
/// Assignment operation
centerline& operator=(const centerline &rhs){
if(this == &rhs) return *this; //test for and handle self-assignment
copy(rhs);
return *this;
|
0d0ef1a9
David Mayerich
Adapted classes t...
|
161
|
}
|
c01958ad
Pavel Govyadinov
merged the two st...
|
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
|
/// 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){
//
// 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 idx; //return the nearest centerline point index
// }
/// Returns the fiber centerline as an array of stim::vec points
std::vector< stim::vec<T> > get_centerline(){
//create an array of stim vectors
std::vector< stim::vec3<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::vec3<T>((T) c[i][0], (T) c[i][1], (T) c[i][2]);
//return the centerline array
return pts;
|
a9ed210f
Pavel Govyadinov
added centerline
|
198
199
200
201
202
|
}
/// Split the fiber at the specified index. If the index is an end point, only one fiber is returned
std::vector< stim::centerline<T> > split(unsigned int idx){
|
c01958ad
Pavel Govyadinov
merged the two st...
|
203
|
std::vector< stim::centerline<T> > fl; //create an array to store up to two fibers
|
a9ed210f
Pavel Govyadinov
added centerline
|
204
205
206
207
208
209
210
211
212
213
|
//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{
|
da89bce9
Jiaming Guo
fixed splitting bugs
|
214
|
unsigned int N1 = idx + 1; //calculate the size of both fibers
|
a9ed210f
Pavel Govyadinov
added centerline
|
215
216
217
218
|
unsigned int N2 = N - idx;
fl.resize(2); //set the array size to 2
|
c01958ad
Pavel Govyadinov
merged the two st...
|
219
220
|
fl[0].init(N1); //set the size of each fiber
fl[1].init(N2);
|
a9ed210f
Pavel Govyadinov
added centerline
|
221
222
|
//copy both halves of the fiber
|
c01958ad
Pavel Govyadinov
merged the two st...
|
223
|
unsigned int i, d;
|
a9ed210f
Pavel Govyadinov
added centerline
|
224
225
|
//first half
|
c01958ad
Pavel Govyadinov
merged the two st...
|
226
227
228
229
|
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
}
|
a9ed210f
Pavel Govyadinov
added centerline
|
230
231
|
//second half
|
c01958ad
Pavel Govyadinov
merged the two st...
|
232
233
234
235
236
|
for(i = 0; i < N2; i++){
for(d = 0; d < 3; d++)
fl[1].c[i][d] = c[idx + i][d];
}
|
a9ed210f
Pavel Govyadinov
added centerline
|
237
238
239
240
241
242
|
}
return fl; //return the array
}
|
c01958ad
Pavel Govyadinov
merged the two st...
|
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
|
/// 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::centerline<T> > connect( stim::centerline<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.
}
}
*/
|
a9ed210f
Pavel Govyadinov
added centerline
|
261
262
263
|
/// Outputs the fiber as a string
std::string str(){
std::stringstream ss;
|
c01958ad
Pavel Govyadinov
merged the two st...
|
264
265
266
267
268
269
270
271
272
|
//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]<<" ";
}
}
|
a9ed210f
Pavel Govyadinov
added centerline
|
273
274
275
|
return ss.str();
}
|
c01958ad
Pavel Govyadinov
merged the two st...
|
276
277
278
|
/// Returns the number of centerline points in the fiber
unsigned int size(){
return N;
|
a9ed210f
Pavel Govyadinov
added centerline
|
279
|
}
|
c72184d1
Jiaming Guo
add many function...
|
280
|
|
a9ed210f
Pavel Govyadinov
added centerline
|
281
|
|
c01958ad
Pavel Govyadinov
merged the two st...
|
282
|
/// Bracket operator returns the element at index i
|
c72184d1
Jiaming Guo
add many function...
|
283
|
|
c01958ad
Pavel Govyadinov
merged the two st...
|
284
285
286
287
|
/// @param i is the index of the element to be returned as a stim::vec
stim::vec<T> operator[](unsigned i){
return get_vec(i);
}
|
c72184d1
Jiaming Guo
add many function...
|
288
|
|
c01958ad
Pavel Govyadinov
merged the two st...
|
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
|
/// Back method returns the last point in the fiber
stim::vec<T> back(){
return get_vec(N-1);
}
////resample a fiber in the network
stim::centerline<T> resample(T spacing)
{
std::cout<<"fiber::resample()"<<std::endl;
std::vector<T> v(3); //v-direction vector of the segment
stim::vec<T> p(3); //- intermediate point to be added
stim::vec<T> p1(3); // p1 - starting point of an segment on the fiber,
stim::vec<T> p2(3); // p2 - ending point,
double sum=0; //distance summation
std::vector<stim::vec<T> > fiberPositions = centerline();
std::vector<stim::vec<T> > newPointList; // initialize list of new resampled points on the fiber
|
a9ed210f
Pavel Govyadinov
added centerline
|
305
|
// for each point on the centerline (skip if it is the last point on centerline)
|
c01958ad
Pavel Govyadinov
merged the two st...
|
306
|
//unsigned int N = fiberPositions.size(); // number of points on the fiber
|
a9ed210f
Pavel Govyadinov
added centerline
|
307
|
for(unsigned int f=0; f< N-1; f++)
|
c01958ad
Pavel Govyadinov
merged the two st...
|
308
|
{
|
a9ed210f
Pavel Govyadinov
added centerline
|
309
|
|
c01958ad
Pavel Govyadinov
merged the two st...
|
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
|
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
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);
}
|
a9ed210f
Pavel Govyadinov
added centerline
|
329
|
}
|
c72184d1
Jiaming Guo
add many function...
|
330
|
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
|
c01958ad
Pavel Govyadinov
merged the two st...
|
331
332
333
334
335
|
newPointList.push_back(fiberPositions[f+1]);
}
newPointList.push_back(fiberPositions[N-1]); //add the last point on the fiber to the new fiber list
centerline newFiber(newPointList);
return newFiber;
|
a9ed210f
Pavel Govyadinov
added centerline
|
336
337
338
339
340
341
342
343
344
345
346
|
}
};
} //end namespace stim
#endif
|