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
|
///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
|
46a9cc26
Jiaming Guo
fix minor errors ...
|
25
|
if (i == size() - 1) return (at(size() - 1) - at(size() - 2)).norm(); //the last direction vector is oriented towards the last line segment
|
a9ed210f
Pavel Govyadinov
added centerline
|
26
|
|
c72184d1
Jiaming Guo
add many function...
|
27
|
//all other direction vectors are the average direction of the two joined line segments
|
83c3121c
David Mayerich
NaN value for |Ax...
|
28
29
30
31
|
vec3<T> a = at(i) - at(i - 1);
vec3<T> b = at(i + 1) - at(i);
vec3<T> ab = a.norm() + b.norm();
return ab.norm();
|
a9ed210f
Pavel Govyadinov
added centerline
|
32
|
}
|
c72184d1
Jiaming Guo
add many function...
|
33
34
35
36
37
38
|
//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
|
39
|
}
|
a9ed210f
Pavel Govyadinov
added centerline
|
40
|
|
c72184d1
Jiaming Guo
add many function...
|
41
42
43
44
|
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
|
45
|
}
|
a9ed210f
Pavel Govyadinov
added centerline
|
46
47
|
}
|
c72184d1
Jiaming Guo
add many function...
|
48
49
50
51
|
void init() {
if (size() == 0) return; //return if there aren't any points
update_L();
}
|
a9ed210f
Pavel Govyadinov
added centerline
|
52
53
54
55
56
|
/// 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...
|
57
58
59
60
61
62
|
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) {
|
0d0ef1a9
David Mayerich
Adapted classes t...
|
63
64
65
66
67
|
for (size_t i = 1; i < L.size(); i++) { //for each point in the centerline
if (L[i] > l) return i - 1; //if we have passed the desired length value, return i-1
}
return L.size() - 1;
/*size_t i = L.size() / 2;
|
c72184d1
Jiaming Guo
add many function...
|
68
69
|
size_t max = L.size() - 1;
size_t min = 0;
|
0d0ef1a9
David Mayerich
Adapted classes t...
|
70
|
while (i < L.size() - 1){
|
c72184d1
Jiaming Guo
add many function...
|
71
72
73
74
75
76
77
78
79
80
81
82
|
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;
}
}
|
0d0ef1a9
David Mayerich
Adapted classes t...
|
83
|
return i;*/
|
c72184d1
Jiaming Guo
add many function...
|
84
|
}
|
a9ed210f
Pavel Govyadinov
added centerline
|
85
|
|
c72184d1
Jiaming Guo
add many function...
|
86
87
88
89
90
91
92
93
94
|
///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
|
95
96
97
98
99
|
}
public:
|
c72184d1
Jiaming Guo
add many function...
|
100
101
|
using std::vector< stim::vec3<T> >::at;
using std::vector< stim::vec3<T> >::size;
|
a9ed210f
Pavel Govyadinov
added centerline
|
102
|
|
c72184d1
Jiaming Guo
add many function...
|
103
104
|
centerline() : std::vector< stim::vec3<T> >() {
init();
|
a9ed210f
Pavel Govyadinov
added centerline
|
105
|
}
|
c72184d1
Jiaming Guo
add many function...
|
106
107
|
centerline(size_t n) : std::vector< stim::vec3<T> >(n){
init();
|
a9ed210f
Pavel Govyadinov
added centerline
|
108
109
|
}
|
c72184d1
Jiaming Guo
add many function...
|
110
111
112
113
|
//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
|
114
115
|
}
|
c72184d1
Jiaming Guo
add many function...
|
116
117
118
119
120
121
|
///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
|
122
|
|
c72184d1
Jiaming Guo
add many function...
|
123
124
125
|
T l = pvalue*L[L.size() - 1];
int idx = findIdx(l);
return p(l, idx);
|
a9ed210f
Pavel Govyadinov
added centerline
|
126
127
|
}
|
0d0ef1a9
David Mayerich
Adapted classes t...
|
128
129
130
131
|
///Update centerline internal parameters (currently the L vector)
void update() {
init();
}
|
c72184d1
Jiaming Guo
add many function...
|
132
133
134
|
///Return the length of the entire centerline
T length() {
return L.back();
|
a9ed210f
Pavel Govyadinov
added centerline
|
135
136
137
138
139
|
}
/// 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...
|
140
141
|
std::vector< stim::centerline<T> > fl; //create an array to store up to two fibers
size_t N = size();
|
a9ed210f
Pavel Govyadinov
added centerline
|
142
143
144
145
146
147
148
149
150
151
|
//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
|
152
|
unsigned int N1 = idx + 1; //calculate the size of both fibers
|
a9ed210f
Pavel Govyadinov
added centerline
|
153
154
155
156
|
unsigned int N2 = N - idx;
fl.resize(2); //set the array size to 2
|
c72184d1
Jiaming Guo
add many function...
|
157
158
|
fl[0] = stim::centerline<T>(N1); //set the size of each fiber
fl[1] = stim::centerline<T>(N2);
|
a9ed210f
Pavel Govyadinov
added centerline
|
159
160
|
//copy both halves of the fiber
|
c72184d1
Jiaming Guo
add many function...
|
161
|
unsigned int i;
|
a9ed210f
Pavel Govyadinov
added centerline
|
162
163
|
//first half
|
c72184d1
Jiaming Guo
add many function...
|
164
165
166
|
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
|
167
168
|
//second half
|
c72184d1
Jiaming Guo
add many function...
|
169
170
171
|
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
|
172
173
174
175
176
177
|
}
return fl; //return the array
}
|
a9ed210f
Pavel Govyadinov
added centerline
|
178
179
180
|
/// Outputs the fiber as a string
std::string str(){
std::stringstream ss;
|
c72184d1
Jiaming Guo
add many function...
|
181
182
183
184
185
|
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
|
186
187
188
|
return ss.str();
}
|
a9ed210f
Pavel Govyadinov
added centerline
|
189
190
|
/// Back method returns the last point in the fiber
|
c72184d1
Jiaming Guo
add many function...
|
191
192
|
stim::vec3<T> back(){
return std::vector< stim::vec3<T> >::back();
|
a9ed210f
Pavel Govyadinov
added centerline
|
193
|
}
|
c72184d1
Jiaming Guo
add many function...
|
194
|
|
a9ed210f
Pavel Govyadinov
added centerline
|
195
196
|
////resample a fiber in the network
stim::centerline<T> resample(T spacing)
|
1599ae65
Jiaming Guo
add swc.h and fun...
|
197
198
|
{
//std::cout<<"fiber::resample()"<<std::endl;
|
a9ed210f
Pavel Govyadinov
added centerline
|
199
|
|
c72184d1
Jiaming Guo
add many function...
|
200
201
202
203
204
205
206
207
208
|
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
|
209
|
// for each point on the centerline (skip if it is the last point on centerline)
|
a9ed210f
Pavel Govyadinov
added centerline
|
210
|
for(unsigned int f=0; f< N-1; f++)
|
c72184d1
Jiaming Guo
add many function...
|
211
212
213
214
|
{
p1 = at(f);
p2 = at(f+1);
v = p2 - p1;
|
a9ed210f
Pavel Govyadinov
added centerline
|
215
|
|
c72184d1
Jiaming Guo
add many function...
|
216
217
218
219
220
221
222
223
224
225
226
|
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
|
227
|
}
|
a9ed210f
Pavel Govyadinov
added centerline
|
228
|
}
|
c72184d1
Jiaming Guo
add many function...
|
229
|
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
|
62d95dcd
Jiaming Guo
fixed resample bugs
|
230
|
new_c.push_back(at(f));
|
c72184d1
Jiaming Guo
add many function...
|
231
232
233
234
|
}
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
|
235
236
237
238
239
240
241
242
243
244
245
|
}
};
} //end namespace stim
#endif
|