5b6c317a
Pavel Govyadinov
implemented the c...
|
1
2
3
4
5
6
7
8
9
10
11
12
13
|
#ifndef STIM_CYLINDER_H
#define STIM_CYLINDER_H
#include <iostream>
#include <stim/math/circle.h>
#include <stim/math/vector.h>
namespace stim
{
template<typename T>
class cylinder
{
private:
|
8495a970
Pavel Govyadinov
added comments to...
|
14
15
16
17
|
stim::circle<T> s; //an arbitrary circle
std::vector< stim::vec<T> > pos; //positions of the cylinder.
std::vector< stim::vec<T> > mags; //radii at each position
std::vector< T > L; //length of the cylinder at each position.
|
9c97e126
David Mayerich
added an axis-ali...
|
18
19
|
///default init
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
20
|
void
|
9c97e126
David Mayerich
added an axis-ali...
|
21
|
init(){
|
10654a1f
Pavel Govyadinov
added the necessa...
|
22
|
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
23
24
|
}
|
8495a970
Pavel Govyadinov
added comments to...
|
25
|
///inits the cylinder from a list of points (inP) and radii (inM)
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
26
|
void
|
9c97e126
David Mayerich
added an axis-ali...
|
27
|
init(std::vector<stim::vec<T> > inP, std::vector<stim::vec<T> > inM){
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
28
29
|
pos = inP;
mags = inM;
|
8495a970
Pavel Govyadinov
added comments to...
|
30
31
|
//calculate each L.
|
10654a1f
Pavel Govyadinov
added the necessa...
|
32
|
L.resize(pos.size()-1);
|
bf23ee36
Pavel Govyadinov
more minor bug fi...
|
33
|
T temp = (T)0;
|
9c97e126
David Mayerich
added an axis-ali...
|
34
|
for(int i = 0; i < L.size(); i++)
|
10654a1f
Pavel Govyadinov
added the necessa...
|
35
|
{
|
bf23ee36
Pavel Govyadinov
more minor bug fi...
|
36
37
|
temp += (pos[i] - pos[i+1]).len();
L[i] = temp;
|
10654a1f
Pavel Govyadinov
added the necessa...
|
38
|
}
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
39
|
}
|
9c97e126
David Mayerich
added an axis-ali...
|
40
|
|
8495a970
Pavel Govyadinov
added comments to...
|
41
|
///returns the direction vector at point idx.
|
10654a1f
Pavel Govyadinov
added the necessa...
|
42
|
stim::vec<T>
|
9c97e126
David Mayerich
added an axis-ali...
|
43
|
d(int idx){
|
10654a1f
Pavel Govyadinov
added the necessa...
|
44
|
return (pos[idx] - pos[idx+1]).norm();
|
10654a1f
Pavel Govyadinov
added the necessa...
|
45
46
|
}
|
8495a970
Pavel Govyadinov
added comments to...
|
47
|
///returns the total length of the line at index j.
|
10654a1f
Pavel Govyadinov
added the necessa...
|
48
|
T
|
9c97e126
David Mayerich
added an axis-ali...
|
49
|
getl(int j){
|
10654a1f
Pavel Govyadinov
added the necessa...
|
50
51
|
for(int i = 0; i < j-1; ++i)
{
|
bf23ee36
Pavel Govyadinov
more minor bug fi...
|
52
53
|
temp += (pos[i] - pos[i+1]).len();
L[i] = temp;
|
10654a1f
Pavel Govyadinov
added the necessa...
|
54
55
56
|
}
}
|
8495a970
Pavel Govyadinov
added comments to...
|
57
58
|
///finds the index of the point closest to the length l on the lower bound.
///binary search.
|
10654a1f
Pavel Govyadinov
added the necessa...
|
59
|
int
|
9c97e126
David Mayerich
added an axis-ali...
|
60
|
findIdx(T l){
|
10654a1f
Pavel Govyadinov
added the necessa...
|
61
|
int i = pos.size()/2;
|
bf23ee36
Pavel Govyadinov
more minor bug fi...
|
62
|
while(i > 0 && i < pos.size())
|
10654a1f
Pavel Govyadinov
added the necessa...
|
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
|
{
if(L[i] < l)
{
i = i/2;
}
else if(L[i] < l && L[i+1] > l)
{
break;
}
else
{
i = i+i/2;
}
}
return i;
}
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
79
|
|
9c97e126
David Mayerich
added an axis-ali...
|
80
81
82
83
84
85
86
87
88
89
90
91
|
//initializes the length array given the current set of positions
void init_length(){
vec<T> p0, p1;
p0 = pos[0]; //initialize the first point in the segment to the first point in the cylinder
T l; //allocate space for the segment length
for(unsigned p = 1; p < pos.size(); p++){ //for each point in the cylinder
p1 = pos[p]; //get the second point in the segment
l = (p1 - p0).len(); //calculate the length of the segment
if(p == 1) L[0] = l; //set the length for the first segment
else L[p-1] = L[p-2] + l; //calculate and set the running length for each additional segment
}
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
92
93
94
|
}
|
9c97e126
David Mayerich
added an axis-ali...
|
95
96
97
98
|
public:
///default constructor
cylinder(){}
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
99
|
///constructor to create a cylinder from a set of points, radii, and the number of sides for the cylinder.
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
100
101
|
///@param inP: Vector of stim vecs composing the points of the centerline.
///@param inM: Vector of stim vecs composing the radii of the centerline.
|
9c97e126
David Mayerich
added an axis-ali...
|
102
|
cylinder(std::vector<stim::vec<T> > inP, std::vector<stim::vec<T> > inM){
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
103
104
105
|
init(inP, inM);
}
|
58ee2b21
David Mayerich
incorporated stim...
|
106
107
108
109
|
///Constructor defines a cylinder with centerline inP and magnitudes of zero
///@param inP: Vector of stim vecs composing the points of the centerline
cylinder(std::vector< stim::vec<T> > inP){
std::vector< stim::vec<T> > inM; //create an array of arbitrary magnitudes
|
9c97e126
David Mayerich
added an axis-ali...
|
110
111
112
113
114
|
stim::vec<T> zero;
zero.push_back(0);
inM.resize(inP.size(), zero); //initialize the magnitude values to zero
|
58ee2b21
David Mayerich
incorporated stim...
|
115
116
117
118
119
120
121
122
123
124
|
init(inP, inM);
}
///Returns the number of points on the cylinder centerline
unsigned int size(){
return pos.size();
}
|
10654a1f
Pavel Govyadinov
added the necessa...
|
125
126
|
///Returns a position vector at the given p-value (p value ranges from 0 to 1).
|
8495a970
Pavel Govyadinov
added comments to...
|
127
128
|
///interpolates the position along the line.
///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
|
10654a1f
Pavel Govyadinov
added the necessa...
|
129
|
stim::vec<T>
|
9c97e126
David Mayerich
added an axis-ali...
|
130
|
p(T pvalue){
|
bf23ee36
Pavel Govyadinov
more minor bug fi...
|
131
132
|
if(pvalue < 0.0 || pvalue > 1.0)
return;
|
10654a1f
Pavel Govyadinov
added the necessa...
|
133
134
|
T l = pvalue*L[L.size()-1];
int idx = findIdx(l);
|
66fe63f3
Pavel Govyadinov
fixed minor typo
|
135
|
return (pos[idx] + (pos[idx+1]-pos[idx])*((l-L[idx])/(L[idx+1]- L[idx])));
|
10654a1f
Pavel Govyadinov
added the necessa...
|
136
137
|
}
|
8495a970
Pavel Govyadinov
added comments to...
|
138
139
140
141
|
///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.
|
10654a1f
Pavel Govyadinov
added the necessa...
|
142
|
stim::vec<T>
|
9c97e126
David Mayerich
added an axis-ali...
|
143
|
p(T l, int idx){
|
66fe63f3
Pavel Govyadinov
fixed minor typo
|
144
|
return (pos[idx] + (pos[idx+1]-pos[idx])*((l-L[idx])/(L[idx+1]- L[idx])));
|
10654a1f
Pavel Govyadinov
added the necessa...
|
145
146
147
|
}
///Returns a radius at the given p-value (p value ranges from 0 to 1).
|
8495a970
Pavel Govyadinov
added comments to...
|
148
149
|
///interpolates the radius along the line.
///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
|
10654a1f
Pavel Govyadinov
added the necessa...
|
150
|
T
|
9c97e126
David Mayerich
added an axis-ali...
|
151
|
r(T pvalue){
|
bf23ee36
Pavel Govyadinov
more minor bug fi...
|
152
153
|
if(pvalue < 0.0 || pvalue > 1.0)
return;
|
10654a1f
Pavel Govyadinov
added the necessa...
|
154
155
|
T l = pvalue*L[L.size()-1];
int idx = findIdx(l);
|
66fe63f3
Pavel Govyadinov
fixed minor typo
|
156
157
|
return (mags[idx] + (mags[idx+1]-mags[idx])*((l-L[idx])/(L[idx+1]- L[idx])));
}
|
10654a1f
Pavel Govyadinov
added the necessa...
|
158
|
|
8495a970
Pavel Govyadinov
added comments to...
|
159
160
161
162
|
///Returns a radius at the given length into the fiber (based on the pvalue).
///Interpolates the position 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.
|
10654a1f
Pavel Govyadinov
added the necessa...
|
163
|
T
|
9c97e126
David Mayerich
added an axis-ali...
|
164
|
r(T l, int idx){
|
66fe63f3
Pavel Govyadinov
fixed minor typo
|
165
166
|
return (mags[idx] + (mags[idx+1]-mags[idx])*((l-L[idx])/(L[idx+1]- L[idx])));
}
|
10654a1f
Pavel Govyadinov
added the necessa...
|
167
|
|
9c97e126
David Mayerich
added an axis-ali...
|
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
|
/// Returns the magnitude at the given index
/// @param i is the index of the desired point
/// @param m is the index of the magnitude value
T ri(unsigned i, unsigned m = 0){
return mags[i][m];
}
/// Adds a new magnitude value to all points
/// @param m is the starting value for the new magnitude
void add_mag(T m = 0){
for(unsigned int p = 0; p < pos.size(); p++)
mags[p].push_back(m);
}
/// Sets a magnitude value
/// @param val is the new value for the magnitude
/// @param p is the point index for the magnitude to be set
/// @param m is the index for the magnitude
void set_mag(T val, unsigned p, unsigned m = 0){
mags[p][m] = val;
}
|
b3a38641
David Mayerich
added the ability...
|
190
191
192
193
194
|
/// Returns the number of magnitude values at each point
unsigned nmags(){
return mags[0].size();
}
|
9c97e126
David Mayerich
added an axis-ali...
|
195
|
|
10654a1f
Pavel Govyadinov
added the necessa...
|
196
197
|
///returns the position of the point with a given pvalue and theta on the surface
|
8495a970
Pavel Govyadinov
added comments to...
|
198
199
200
|
///in x, y, z coordinates. Theta is in degrees from 0 to 360.
///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
///@param theta: the angle to the point of a circle.
|
10654a1f
Pavel Govyadinov
added the necessa...
|
201
202
203
|
stim::vec<T>
surf(T pvalue, T theta)
{
|
bf23ee36
Pavel Govyadinov
more minor bug fi...
|
204
205
|
if(pvalue < 0.0 || pvalue > 1.0)
return;
|
10654a1f
Pavel Govyadinov
added the necessa...
|
206
207
|
T l = pvalue*L[L.size()-1];
int idx = findIdx(l);
|
9c97e126
David Mayerich
added an axis-ali...
|
208
|
stim::vec<T> ps = p(l, idx);
|
10654a1f
Pavel Govyadinov
added the necessa...
|
209
210
211
212
213
214
|
T m = r(l, idx);
stim::vec<T> dr = d(idx);
s = stim::circle<T>(ps, m, dr);
return(s.p(theta));
}
|
8495a970
Pavel Govyadinov
added comments to...
|
215
|
///returns a vector of points necessary to create a circle at every position in the fiber.
|
9c97e126
David Mayerich
added an axis-ali...
|
216
|
///@param sides: the number of sides of each circle.
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
|
std::vector<std::vector<vec<T> > >
getPoints(int sides)
{
if(pos.size() < 2)
{
return;
} else {
std::vector<std::vector <vec<T> > > points;
points.resize(pos.size());
stim::vec<T> d = (pos[0] - pos[1]).norm();
s = stim::circle<T>(pos[0], mags[0][0], d);
points[0] = s.getPoints(sides);
for(int i = 1; i < pos.size(); i++)
{
d = (pos[i] - pos[i-1]).norm();
s.center(pos[i]);
s.normal(d);
s.scale(mags[i][0]/mags[i-1][0], mags[i][0]/mags[i-1][0]);
points[i] = s.getPoints(sides);
}
return points;
}
}
|
58ee2b21
David Mayerich
incorporated stim...
|
240
241
242
243
|
/// Allows a point on the centerline to be accessed using bracket notation
vec<T> operator[](unsigned int i){
|
58ee2b21
David Mayerich
incorporated stim...
|
244
|
return pos[i];
|
58ee2b21
David Mayerich
incorporated stim...
|
245
246
247
248
249
250
251
|
}
/// Returns the total length of the cylinder centerline
T length(){
return L.back();
}
|
9c97e126
David Mayerich
added an axis-ali...
|
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
|
/// Integrates a magnitude value along the cylinder.
/// @param m is the magnitude value to be integrated (this is usually the radius)
T integrate(unsigned m = 0){
T M = 0; //initialize the integral to zero
T m0, m1; //allocate space for both magnitudes in a single segment
//vec<T> p0, p1; //allocate space for both points in a single segment
m0 = mags[0][m]; //initialize the first point and magnitude to the first point in the cylinder
//p0 = pos[0];
T len = L[0]; //allocate space for the segment length
//for every consecutive point in the cylinder
for(unsigned p = 1; p < pos.size(); p++){
//p1 = pos[p]; //get the position and magnitude for the next point
m1 = mags[p][m];
if(p > 1) len = (L[p-1] - L[p-2]); //calculate the segment length using the L array
//add the average magnitude, weighted by the segment length
M += (m0 + m1)/2.0 * len;
m0 = m1; //move to the next segment by shifting points
}
return M; //return the integral
}
/// Averages a magnitude value across the cylinder
/// @param m is the magnitude value to be averaged (this is usually the radius)
T average(unsigned m = 0){
//return the average magnitude
return integrate(m) / L.back();
}
|
58ee2b21
David Mayerich
incorporated stim...
|
290
291
292
293
294
295
|
/// Resamples the cylinder to provide a maximum distance of "spacing" between centerline points. All current
/// centerline points are guaranteed to exist in the new cylinder
/// @param spacing is the maximum spacing allowed between sample points
cylinder<T> resample(T spacing){
std::vector< vec<T> > result;
|
9c97e126
David Mayerich
added an axis-ali...
|
296
|
|
58ee2b21
David Mayerich
incorporated stim...
|
297
298
299
300
301
302
303
304
305
306
307
|
vec<T> p0 = pos[0]; //initialize p0 to the first point on the centerline
vec<T> p1;
unsigned N = size(); //number of points in the current centerline
//for each line segment on the centerline
for(unsigned int i = 1; i < N; i++){
p1 = pos[i]; //get the second point in the line segment
vec<T> v = p1 - p0; //calculate the vector between these two points
T d = v.len(); //calculate the distance between these two points (length of the line segment)
|
9c97e126
David Mayerich
added an axis-ali...
|
308
|
unsigned nsteps = d / spacing+1; //calculate the number of steps to take along the segment to meet the spacing criteria
|
58ee2b21
David Mayerich
incorporated stim...
|
309
310
311
312
313
314
315
316
317
318
319
320
|
T stepsize = 1.0 / nsteps; //calculate the parametric step size between new centerline points
//for each step along the line segment
for(unsigned s = 0; s < nsteps; s++){
T alpha = stepsize * s; //calculate the fraction of the distance along the line segment covered
result.push_back(p0 + alpha * v); //push the point at alpha position along the line segment
}
p0 = p1; //shift the points to move to the next line segment
}
result.push_back(pos[size() - 1]); //push the last point in the centerline
|
9c97e126
David Mayerich
added an axis-ali...
|
321
|
|
58ee2b21
David Mayerich
incorporated stim...
|
322
323
324
|
return cylinder<T>(result);
}
|
9c97e126
David Mayerich
added an axis-ali...
|
325
|
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
326
327
328
329
|
};
}
#endif
|