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>
|
a9ed210f
Pavel Govyadinov
added centerline
|
6
7
8
9
10
11
12
13
|
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>
|
c72184d1
Jiaming Guo
add many function...
|
14
|
class centerline : public std::vector< stim::vec3<T> >{
|
a9ed210f
Pavel Govyadinov
added centerline
|
15
16
|
protected:
|
a9ed210f
Pavel Govyadinov
added centerline
|
17
|
|
c72184d1
Jiaming Guo
add many function...
|
18
|
std::vector<T> L; //stores the integrated length along the fiber (used for parameterization)
|
a9ed210f
Pavel Govyadinov
added centerline
|
19
|
|
c72184d1
Jiaming Guo
add many function...
|
20
21
22
23
24
25
|
///Return the normalized direction vector at point i (average of the incoming and outgoing directions)
vec3<T> d(size_t i) {
if (size() <= 1) return vec3<T>(0, 0, 0); //if there is insufficient information to calculate the direction, return a null vector
if (size() == 2) return (at(1) - at(0)).norm(); //if there are only two points, the direction vector at both is the direction of the line segment
if (i == 0) return (at(1) - at(0)).norm(); //the first direction vector is oriented towards the first line segment
if (i == size() - 1) return at(size() - 1) - at(size() - 2); //the last direction vector is oriented towards the last line segment
|
a9ed210f
Pavel Govyadinov
added centerline
|
26
|
|
c72184d1
Jiaming Guo
add many function...
|
27
28
|
//all other direction vectors are the average direction of the two joined line segments
return ((at(i) - at(i - 1)).norm() + (at(i + 1) - at(i)).norm()).norm();
|
a9ed210f
Pavel Govyadinov
added centerline
|
29
|
}
|
c72184d1
Jiaming Guo
add many function...
|
30
31
32
33
34
35
|
//initializes the integrated length vector to make parameterization easier, starting with index idx (all previous indices are assumed to be correct)
void update_L(size_t start = 0) {
L.resize(size()); //allocate space for the L array
if (start == 0) {
L[0] = 0; //initialize the length value for the first point to zero (0)
start++;
|
a9ed210f
Pavel Govyadinov
added centerline
|
36
|
}
|
a9ed210f
Pavel Govyadinov
added centerline
|
37
|
|
c72184d1
Jiaming Guo
add many function...
|
38
39
40
41
|
stim::vec3<T> d;
for (size_t i = start; i < size(); i++) { //for each line segment in the centerline
d = at(i) - at(i - 1);
L[i] = L[i - 1] + d.len(); //calculate the running length total
|
a9ed210f
Pavel Govyadinov
added centerline
|
42
|
}
|
a9ed210f
Pavel Govyadinov
added centerline
|
43
44
|
}
|
c72184d1
Jiaming Guo
add many function...
|
45
46
47
48
|
void init() {
if (size() == 0) return; //return if there aren't any points
update_L();
}
|
a9ed210f
Pavel Govyadinov
added centerline
|
49
50
51
52
53
|
/// 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){
|
c72184d1
Jiaming Guo
add many function...
|
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
|
return std::vector< stim::vec3<T> >::at(i);
}
///finds the index of the point closest to the length l on the lower bound.
///binary search.
size_t findIdx(T l) {
size_t i = L.size() / 2;
size_t max = L.size() - 1;
size_t min = 0;
while (i > 0 && i < L.size() - 1){
if (l < L[i]) {
max = i;
i = min + (max - min) / 2;
}
else if (L[i] <= l && L[i + 1] >= l) {
break;
}
else {
min = i;
i = min + (max - min) / 2;
}
}
return i;
}
|
a9ed210f
Pavel Govyadinov
added centerline
|
78
|
|
c72184d1
Jiaming Guo
add many function...
|
79
80
81
82
83
84
85
86
87
|
///Returns a position vector at the given length into the fiber (based on the pvalue).
///Interpolates the radius along the line.
///@param l: the location of the in the cylinder.
///@param idx: integer location of the point closest to l but prior to it.
stim::vec3<T> p(T l, int idx) {
T rat = (l - L[idx]) / (L[idx + 1] - L[idx]);
stim::vec3<T> v1 = at(idx);
stim::vec3<T> v2 = at(idx + 1);
return(v1 + (v2 - v1)*rat);
|
a9ed210f
Pavel Govyadinov
added centerline
|
88
89
90
91
92
|
}
public:
|
c72184d1
Jiaming Guo
add many function...
|
93
94
|
using std::vector< stim::vec3<T> >::at;
using std::vector< stim::vec3<T> >::size;
|
a9ed210f
Pavel Govyadinov
added centerline
|
95
|
|
c72184d1
Jiaming Guo
add many function...
|
96
97
|
centerline() : std::vector< stim::vec3<T> >() {
init();
|
a9ed210f
Pavel Govyadinov
added centerline
|
98
|
}
|
c72184d1
Jiaming Guo
add many function...
|
99
100
|
centerline(size_t n) : std::vector< stim::vec3<T> >(n){
init();
|
a9ed210f
Pavel Govyadinov
added centerline
|
101
102
|
}
|
c72184d1
Jiaming Guo
add many function...
|
103
104
105
106
|
//overload the push_back function to update the length vector
void push_back(stim::vec3<T> p) {
std::vector< stim::vec3<T> >::push_back(p);
update_L(size() - 1);
|
a9ed210f
Pavel Govyadinov
added centerline
|
107
108
|
}
|
c72184d1
Jiaming Guo
add many function...
|
109
110
111
112
113
114
|
///Returns a position vector at the given p-value (p value ranges from 0 to 1).
///interpolates the position along the line.
///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
stim::vec3<T> p(T pvalue) {
if (pvalue <= 0.0) return at(0); //return the first element
if (pvalue >= 1.0) return back(); //return the last element
|
a9ed210f
Pavel Govyadinov
added centerline
|
115
|
|
c72184d1
Jiaming Guo
add many function...
|
116
117
118
|
T l = pvalue*L[L.size() - 1];
int idx = findIdx(l);
return p(l, idx);
|
a9ed210f
Pavel Govyadinov
added centerline
|
119
120
|
}
|
c72184d1
Jiaming Guo
add many function...
|
121
122
123
|
///Return the length of the entire centerline
T length() {
return L.back();
|
a9ed210f
Pavel Govyadinov
added centerline
|
124
125
126
127
128
|
}
/// 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){
|
c72184d1
Jiaming Guo
add many function...
|
129
130
|
std::vector< stim::centerline<T> > fl; //create an array to store up to two fibers
size_t N = size();
|
a9ed210f
Pavel Govyadinov
added centerline
|
131
132
133
134
135
136
137
138
139
140
|
//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{
|
c72184d1
Jiaming Guo
add many function...
|
141
|
unsigned int N1 = idx; //calculate the size of both fibers
|
a9ed210f
Pavel Govyadinov
added centerline
|
142
143
144
145
|
unsigned int N2 = N - idx;
fl.resize(2); //set the array size to 2
|
c72184d1
Jiaming Guo
add many function...
|
146
147
|
fl[0] = stim::centerline<T>(N1); //set the size of each fiber
fl[1] = stim::centerline<T>(N2);
|
a9ed210f
Pavel Govyadinov
added centerline
|
148
149
|
//copy both halves of the fiber
|
c72184d1
Jiaming Guo
add many function...
|
150
|
unsigned int i;
|
a9ed210f
Pavel Govyadinov
added centerline
|
151
152
|
//first half
|
c72184d1
Jiaming Guo
add many function...
|
153
154
155
|
for(i = 0; i < N1; i++) //for each centerline point
fl[0][i] = std::vector< stim::vec3<T> >::at(i);
fl[0].init(); //initialize the length vector
|
a9ed210f
Pavel Govyadinov
added centerline
|
156
157
|
//second half
|
c72184d1
Jiaming Guo
add many function...
|
158
159
160
|
for(i = 0; i < N2; i++)
fl[1][i] = std::vector< stim::vec3<T> >::at(idx+i);
fl[1].init(); //initialize the length vector
|
a9ed210f
Pavel Govyadinov
added centerline
|
161
162
163
164
165
166
|
}
return fl; //return the array
}
|
a9ed210f
Pavel Govyadinov
added centerline
|
167
168
169
|
/// Outputs the fiber as a string
std::string str(){
std::stringstream ss;
|
c72184d1
Jiaming Guo
add many function...
|
170
171
172
173
174
|
size_t N = std::vector< stim::vec3<T> >::size();
ss << "---------[" << N << "]---------" << std::endl;
for (size_t i = 0; i < N; i++)
ss << std::vector< stim::vec3<T> >::at(i) << std::endl;
ss << "--------------------" << std::endl;
|
a9ed210f
Pavel Govyadinov
added centerline
|
175
176
177
|
return ss.str();
}
|
a9ed210f
Pavel Govyadinov
added centerline
|
178
179
|
/// Back method returns the last point in the fiber
|
c72184d1
Jiaming Guo
add many function...
|
180
181
|
stim::vec3<T> back(){
return std::vector< stim::vec3<T> >::back();
|
a9ed210f
Pavel Govyadinov
added centerline
|
182
|
}
|
c72184d1
Jiaming Guo
add many function...
|
183
|
|
a9ed210f
Pavel Govyadinov
added centerline
|
184
185
186
187
188
|
////resample a fiber in the network
stim::centerline<T> resample(T spacing)
{
std::cout<<"fiber::resample()"<<std::endl;
|
c72184d1
Jiaming Guo
add many function...
|
189
190
191
192
193
194
195
196
197
|
stim::vec3<T> v; //v-direction vector of the segment
stim::vec3<T> p; //- intermediate point to be added
stim::vec3<T> p1; // p1 - starting point of an segment on the fiber,
stim::vec3<T> p2; // p2 - ending point,
//double sum=0; //distance summation
size_t N = size();
centerline<T> new_c; // initialize list of new resampled points on the fiber
|
a9ed210f
Pavel Govyadinov
added centerline
|
198
|
// for each point on the centerline (skip if it is the last point on centerline)
|
a9ed210f
Pavel Govyadinov
added centerline
|
199
|
for(unsigned int f=0; f< N-1; f++)
|
c72184d1
Jiaming Guo
add many function...
|
200
201
202
203
|
{
p1 = at(f);
p2 = at(f+1);
v = p2 - p1;
|
a9ed210f
Pavel Govyadinov
added centerline
|
204
|
|
c72184d1
Jiaming Guo
add many function...
|
205
206
207
208
209
210
211
212
213
214
215
|
T lengthSegment = v.len(); //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
p = p1 + v * (step / lengthSegment);
// add this resampled points to the new fiber list
new_c.push_back(p);
|
a9ed210f
Pavel Govyadinov
added centerline
|
216
|
}
|
a9ed210f
Pavel Govyadinov
added centerline
|
217
|
}
|
c72184d1
Jiaming Guo
add many function...
|
218
219
220
221
222
223
|
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
new_c.push_back(at(f+1));
}
new_c.push_back(at(N-1)); //add the last point on the fiber to the new fiber list
//centerline newFiber(newPointList);
return new_c;
|
a9ed210f
Pavel Govyadinov
added centerline
|
224
225
226
227
228
229
230
231
232
233
234
|
}
};
} //end namespace stim
#endif
|