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
|
stim::circle<T> s; //an arbitrary circle
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
15
|
std::vector<stim::circle<T> > e;
|
3e5d3ad3
Pavel Govyadinov
merged the change...
|
16
|
std::vector<stim::vec<T> > mags;
|
8495a970
Pavel Govyadinov
added comments to...
|
17
|
std::vector< T > L; //length of the cylinder at each position.
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
18
19
|
///default init
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
20
|
void
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
21
22
|
init()
{
|
10654a1f
Pavel Govyadinov
added the necessa...
|
23
|
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
24
25
|
}
|
8495a970
Pavel Govyadinov
added comments to...
|
26
|
///inits the cylinder from a list of points (inP) and radii (inM)
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
27
|
void
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
28
29
|
init(std::vector<stim::vec<T> > inP, std::vector<stim::vec<T> > inM)
{
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
30
|
mags = inM;
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
31
32
33
34
35
|
stim::vec<float> v1;
stim::vec<float> v2;
e.resize(inP.size());
if(inP.size() < 2)
return;
|
8495a970
Pavel Govyadinov
added comments to...
|
36
37
|
//calculate each L.
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
38
|
L.resize(inP.size());
|
bf23ee36
Pavel Govyadinov
more minor bug fi...
|
39
|
T temp = (T)0;
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
40
41
|
L[0] = 0;
for(int i = 1; i < L.size(); i++)
|
10654a1f
Pavel Govyadinov
added the necessa...
|
42
|
{
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
43
|
temp += (inP[i-1] - inP[i]).len();
|
bf23ee36
Pavel Govyadinov
more minor bug fi...
|
44
|
L[i] = temp;
|
10654a1f
Pavel Govyadinov
added the necessa...
|
45
|
}
|
9c97e126
David Mayerich
added an axis-ali...
|
46
|
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
|
stim::vec<T> dr = (inP[1] - inP[0]).norm();
s = stim::circle<T>(inP[0], inM[0][0], dr, stim::vec<T>(1,0,0));
e[0] = s;
for(int i = 1; i < inP.size()-1; i++)
{
s.center(inP[i]);
v1 = (inP[i] - inP[i-1]).norm();
v2 = (inP[i+1] - inP[i]).norm();
dr = (v1+v2).norm();
s.normal(dr);
s.scale(inM[i][0]/inM[i-1][0]);
e[i] = s;
}
int j = inP.size()-1;
s.center(inP[j]);
dr = (inP[j] - inP[j-1]).norm();
s.normal(dr);
s.scale(inM[j][0]/inM[j-1][0]);
e[j] = s;
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
67
68
|
}
|
8495a970
Pavel Govyadinov
added comments to...
|
69
|
///returns the direction vector at point idx.
|
10654a1f
Pavel Govyadinov
added the necessa...
|
70
|
stim::vec<T>
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
71
72
73
74
|
d(int idx)
{
if(idx == 0)
{
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
75
|
return (e[idx+1].P - e[idx].P).norm();
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
76
|
}
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
77
|
else if(idx == e.size()-1)
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
78
|
{
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
79
|
return (e[idx].P - e[idx-1].P).norm();
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
80
81
82
|
}
else
{
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
83
84
85
|
// return (e[idx+1].P - e[idx].P).norm();
stim::vec<float> v1 = (e[idx].P-e[idx-1].P).norm();
stim::vec<float> v2 = (e[idx+1].P-e[idx].P).norm();
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
86
|
return (v1+v2).norm();
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
87
88
89
|
}
// return e[idx].N;
|
10654a1f
Pavel Govyadinov
added the necessa...
|
90
91
|
}
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
92
93
|
stim::vec<T>
d(T l, int idx)
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
94
|
{
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
95
|
if(idx == 0 || idx == e.size()-1)
|
10654a1f
Pavel Govyadinov
added the necessa...
|
96
|
{
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
97
|
return e[idx].N;
|
10654a1f
Pavel Govyadinov
added the necessa...
|
98
|
}
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
99
100
101
102
103
|
else
{
T rat = (l-L[idx])/(L[idx+1]-L[idx]);
return( e[idx].N + (e[idx+1].N - e[idx].N)*rat);
}
|
10654a1f
Pavel Govyadinov
added the necessa...
|
104
105
|
}
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
106
|
|
8495a970
Pavel Govyadinov
added comments to...
|
107
108
|
///finds the index of the point closest to the length l on the lower bound.
///binary search.
|
10654a1f
Pavel Govyadinov
added the necessa...
|
109
|
int
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
110
111
112
113
114
|
findIdx(T l)
{
unsigned int i = L.size()/2;
unsigned int max = L.size()-1;
unsigned int min = 0;
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
115
|
while(i > 0 && i < L.size()-1)
|
10654a1f
Pavel Govyadinov
added the necessa...
|
116
|
{
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
117
118
|
// std::cerr << "Trying " << i << std::endl;
// std::cerr << "l is " << l << ", L[" << i << "]" << L[i] << std::endl;
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
119
|
if(l < L[i])
|
10654a1f
Pavel Govyadinov
added the necessa...
|
120
|
{
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
121
|
max = i;
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
122
|
i = min+(max-min)/2;
|
10654a1f
Pavel Govyadinov
added the necessa...
|
123
|
}
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
124
|
else if(L[i] <= l && L[i+1] >= l)
|
10654a1f
Pavel Govyadinov
added the necessa...
|
125
126
127
128
129
|
{
break;
}
else
{
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
130
|
min = i;
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
131
|
i = min+(max-min)/2;
|
10654a1f
Pavel Govyadinov
added the necessa...
|
132
133
134
135
|
}
}
return i;
}
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
136
|
|
9c97e126
David Mayerich
added an axis-ali...
|
137
138
|
public:
///default constructor
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
139
140
141
142
|
cylinder()
{
}
|
9c97e126
David Mayerich
added an axis-ali...
|
143
|
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
144
|
///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...
|
145
146
|
///@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...
|
147
|
cylinder(std::vector<stim::vec<T> > inP, std::vector<stim::vec<T> > inM){
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
148
149
150
|
init(inP, inM);
}
|
58ee2b21
David Mayerich
incorporated stim...
|
151
152
153
154
|
///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...
|
155
156
157
158
159
|
stim::vec<T> zero;
zero.push_back(0);
inM.resize(inP.size(), zero); //initialize the magnitude values to zero
|
58ee2b21
David Mayerich
incorporated stim...
|
160
161
162
163
164
165
166
167
168
169
|
init(inP, inM);
}
///Returns the number of points on the cylinder centerline
unsigned int size(){
return pos.size();
}
|
10654a1f
Pavel Govyadinov
added the necessa...
|
170
171
|
///Returns a position vector at the given p-value (p value ranges from 0 to 1).
|
8495a970
Pavel Govyadinov
added comments to...
|
172
173
|
///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...
|
174
|
stim::vec<T>
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
175
176
|
p(T pvalue)
{
|
bf23ee36
Pavel Govyadinov
more minor bug fi...
|
177
|
if(pvalue < 0.0 || pvalue > 1.0)
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
178
179
180
|
{
return stim::vec<float>(-1,-1,-1);
}
|
10654a1f
Pavel Govyadinov
added the necessa...
|
181
182
|
T l = pvalue*L[L.size()-1];
int idx = findIdx(l);
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
183
|
T rat = (l-L[idx])/(L[idx+1]-L[idx]);
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
184
|
return( e[idx].P + (e[idx+1].P-e[idx].P)*rat);
|
10654a1f
Pavel Govyadinov
added the necessa...
|
185
186
|
}
|
8495a970
Pavel Govyadinov
added comments to...
|
187
188
189
190
|
///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...
|
191
|
stim::vec<T>
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
192
193
194
|
p(T l, int idx)
{
T rat = (l-L[idx])/(L[idx+1]-L[idx]);
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
195
|
return( e[idx].P + (e[idx+1].P-e[idx].P)*rat);
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
196
197
|
// return(
// return (pos[idx] + (pos[idx+1]-pos[idx])*((l-L[idx])/(L[idx+1]- L[idx])));
|
10654a1f
Pavel Govyadinov
added the necessa...
|
198
199
200
|
}
///Returns a radius at the given p-value (p value ranges from 0 to 1).
|
8495a970
Pavel Govyadinov
added comments to...
|
201
202
|
///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...
|
203
|
T
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
204
205
|
r(T pvalue)
{
|
bf23ee36
Pavel Govyadinov
more minor bug fi...
|
206
207
|
if(pvalue < 0.0 || pvalue > 1.0)
return;
|
10654a1f
Pavel Govyadinov
added the necessa...
|
208
209
|
T l = pvalue*L[L.size()-1];
int idx = findIdx(l);
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
210
|
return (e[idx].U.len() + (e[idx+1].U.len() - e[idx].U.len())*((l-L[idx])/(L[idx+1]- L[idx])));
|
66fe63f3
Pavel Govyadinov
fixed minor typo
|
211
|
}
|
10654a1f
Pavel Govyadinov
added the necessa...
|
212
|
|
8495a970
Pavel Govyadinov
added comments to...
|
213
214
215
216
|
///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...
|
217
|
T
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
218
219
220
|
r(T l, int idx)
{
T rat = (l-L[idx])/(L[idx+1]-L[idx]);
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
221
|
return( e[idx].U.len() + (e[idx+1].U.len() - e[idx].U.len())*rat);
|
66fe63f3
Pavel Govyadinov
fixed minor typo
|
222
|
}
|
10654a1f
Pavel Govyadinov
added the necessa...
|
223
|
|
9c97e126
David Mayerich
added an axis-ali...
|
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
|
/// 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...
|
246
247
248
249
250
|
/// Returns the number of magnitude values at each point
unsigned nmags(){
return mags[0].size();
}
|
10654a1f
Pavel Govyadinov
added the necessa...
|
251
|
///returns the position of the point with a given pvalue and theta on the surface
|
8495a970
Pavel Govyadinov
added comments to...
|
252
253
254
|
///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...
|
255
256
257
|
stim::vec<T>
surf(T pvalue, T theta)
{
|
bf23ee36
Pavel Govyadinov
more minor bug fi...
|
258
|
if(pvalue < 0.0 || pvalue > 1.0)
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
259
260
261
|
{
return stim::vec<float>(-1,-1,-1);
} else {
|
10654a1f
Pavel Govyadinov
added the necessa...
|
262
263
|
T l = pvalue*L[L.size()-1];
int idx = findIdx(l);
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
264
|
stim::vec<T> ps = p(l, idx);
|
10654a1f
Pavel Govyadinov
added the necessa...
|
265
|
T m = r(l, idx);
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
266
267
268
269
|
s = e[idx];
s.center(ps);
s.normal(d(l, idx));
s.scale(m/e[idx].U.len());
|
10654a1f
Pavel Govyadinov
added the necessa...
|
270
|
return(s.p(theta));
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
271
|
}
|
10654a1f
Pavel Govyadinov
added the necessa...
|
272
273
|
}
|
8495a970
Pavel Govyadinov
added comments to...
|
274
|
///returns a vector of points necessary to create a circle at every position in the fiber.
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
275
|
///@param sides: the number of sides of each circle.
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
276
277
278
|
std::vector<std::vector<vec<T> > >
getPoints(int sides)
{
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
279
280
281
|
std::vector<std::vector <vec<T> > > points;
points.resize(e.size());
for(int i = 0; i < e.size(); i++)
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
282
|
{
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
283
|
points[i] = e[i].getPoints(sides);
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
284
|
}
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
285
|
return points;
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
286
|
}
|
58ee2b21
David Mayerich
incorporated stim...
|
287
|
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
288
289
290
291
292
293
|
///returns the total length of the line at index j.
T
getl(int j)
{
return (L[j]);
}
|
58ee2b21
David Mayerich
incorporated stim...
|
294
295
296
|
/// Allows a point on the centerline to be accessed using bracket notation
vec<T> operator[](unsigned int i){
|
3e5d3ad3
Pavel Govyadinov
merged the change...
|
297
|
return e[i].P;
|
58ee2b21
David Mayerich
incorporated stim...
|
298
299
300
301
302
303
304
|
}
/// Returns the total length of the cylinder centerline
T length(){
return L.back();
}
|
9c97e126
David Mayerich
added an axis-ali...
|
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
|
/// 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...
|
343
344
345
346
347
348
|
/// 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...
|
349
|
|
3e5d3ad3
Pavel Govyadinov
merged the change...
|
350
|
vec<T> p0 = e[i].P; //initialize p0 to the first point on the centerline
|
58ee2b21
David Mayerich
incorporated stim...
|
351
352
353
354
355
|
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++){
|
3e5d3ad3
Pavel Govyadinov
merged the change...
|
356
|
p1 = e[i].P; //get the second point in the line segment
|
58ee2b21
David Mayerich
incorporated stim...
|
357
358
359
360
|
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...
|
361
|
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...
|
362
363
364
365
366
367
368
369
370
371
372
|
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
}
|
3e5d3ad3
Pavel Govyadinov
merged the change...
|
373
|
result.push_back(e[size() - 1].P); //push the last point in the centerline
|
9c97e126
David Mayerich
added an axis-ali...
|
374
|
|
58ee2b21
David Mayerich
incorporated stim...
|
375
376
377
|
return cylinder<T>(result);
}
|
9c97e126
David Mayerich
added an axis-ali...
|
378
|
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
379
|
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
380
381
382
383
|
};
}
#endif
|