5b6c317a
Pavel Govyadinov
implemented the c...
|
1
2
3
4
|
#ifndef STIM_CYLINDER_H
#define STIM_CYLINDER_H
#include <iostream>
#include <stim/math/circle.h>
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
5
|
#include <stim/biomodels/centerline.h>
|
e21d1051
David Mayerich
cylinder can buil...
|
6
|
#include <stim/visualization/obj.h>
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
7
|
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
8
9
10
11
|
namespace stim
{
template<typename T>
|
e21d1051
David Mayerich
cylinder can buil...
|
12
13
|
class cylinder : public centerline<T> {
protected:
|
ee17b90b
Jiaming Guo
make it compatibl...
|
14
15
|
using stim::centerline<T>::d;
|
e21d1051
David Mayerich
cylinder can buil...
|
16
|
|
e21d1051
David Mayerich
cylinder can buil...
|
17
|
std::vector< stim::vec3<T> > U; //stores the array of U vectors defining the Frenet frame
|
7115e973
David Mayerich
removed the magni...
|
18
|
std::vector< T > R; //stores a list of magnitudes for each point in the centerline (assuming mags[0] is the radius)
|
e21d1051
David Mayerich
cylinder can buil...
|
19
|
|
e21d1051
David Mayerich
cylinder can buil...
|
20
|
using stim::centerline<T>::findIdx;
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
21
|
|
e21d1051
David Mayerich
cylinder can buil...
|
22
23
24
25
|
//calculates the U values for each point to initialize the frenet frame
// this function assumes that the centerline has already been set
void init() {
U.resize(size()); //allocate space for the frenet frame vectors
|
0d0ef1a9
David Mayerich
Adapted classes t...
|
26
|
R.resize(size());
|
e21d1051
David Mayerich
cylinder can buil...
|
27
28
29
30
31
32
|
stim::circle<T> c; //create a circle
stim::vec3<T> d0, d1;
for (size_t i = 0; i < size() - 1; i++) { //for each line segment in the centerline
c.rotate(d(i)); //rotate the circle to match that normal
U[i] = c.U; //save the U vector from the circle
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
33
|
}
|
e21d1051
David Mayerich
cylinder can buil...
|
34
35
36
|
U[size() - 1] = c.U; //for the last point, duplicate the final frenet frame vector
}
|
e21d1051
David Mayerich
cylinder can buil...
|
37
|
public:
|
e21d1051
David Mayerich
cylinder can buil...
|
38
|
|
eb554b48
David Mayerich
bug fixes for Net...
|
39
40
|
using stim::centerline<T>::size;
using stim::centerline<T>::at;
|
ee17b90b
Jiaming Guo
make it compatibl...
|
41
42
|
using stim::centerline<T>::L;
using stim::centerline<T>::length;
|
eb554b48
David Mayerich
bug fixes for Net...
|
43
|
|
ee17b90b
Jiaming Guo
make it compatibl...
|
44
|
cylinder() : centerline<T>(){}
|
e21d1051
David Mayerich
cylinder can buil...
|
45
|
|
ee17b90b
Jiaming Guo
make it compatibl...
|
46
|
cylinder(centerline<T>c) : centerline<T>(c) {
|
e21d1051
David Mayerich
cylinder can buil...
|
47
|
init();
|
e21d1051
David Mayerich
cylinder can buil...
|
48
49
|
}
|
ee17b90b
Jiaming Guo
make it compatibl...
|
50
51
52
53
54
|
//cylinder(centerline<T>c, T r) : centerline(c) {
// init();
// //add_mag(r);
//}
|
e21d1051
David Mayerich
cylinder can buil...
|
55
|
//initialize a cylinder with a list of points and magnitude values
|
ee17b90b
Jiaming Guo
make it compatibl...
|
56
57
58
59
|
//cylinder(centerline<T>c, std::vector<T> r) : centerline(c){
// init();
// //add_mag(r);
//}
|
e21d1051
David Mayerich
cylinder can buil...
|
60
|
|
58674036
Jiaming Guo
fix minor error f...
|
61
62
63
64
65
66
|
//copy the original radius
void copy_r(std::vector<T> radius) {
for (unsigned i = 0; i < radius.size(); i++)
R[i] = radius[i];
}
|
e21d1051
David Mayerich
cylinder can buil...
|
67
68
69
70
|
///Returns magnitude i 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.
|
7115e973
David Mayerich
removed the magni...
|
71
|
T r(T l, int idx) {
|
e21d1051
David Mayerich
cylinder can buil...
|
72
|
T a = (l - L[idx]) / (L[idx + 1] - L[idx]);
|
7115e973
David Mayerich
removed the magni...
|
73
|
T v2 = (R[idx] + (R[idx + 1] - R[idx])*a);
|
e21d1051
David Mayerich
cylinder can buil...
|
74
75
76
77
78
79
80
|
return v2;
}
///Returns the ith magnitude at the given p-value (p value ranges from 0 to 1).
///interpolates the radius along the line.
///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
|
720240c8
Jiaming Guo
everything now is...
|
81
|
T rl(T pvalue) {
|
7115e973
David Mayerich
removed the magni...
|
82
83
|
if (pvalue <= 0.0) return R[0];
if (pvalue >= 1.0) return R[size() - 1];
|
e21d1051
David Mayerich
cylinder can buil...
|
84
85
86
|
T l = pvalue*L[L.size() - 1];
int idx = findIdx(l);
|
0d0ef1a9
David Mayerich
Adapted classes t...
|
87
|
return r(l, idx);
|
e21d1051
David Mayerich
cylinder can buil...
|
88
89
|
}
|
eb554b48
David Mayerich
bug fixes for Net...
|
90
91
|
/// Returns the magnitude at the given index
/// @param i is the index of the desired point
|
7115e973
David Mayerich
removed the magni...
|
92
|
/// @param r is the index of the magnitude value
|
720240c8
Jiaming Guo
everything now is...
|
93
94
95
|
T r(unsigned i) {
return R[i];
}
|
eb554b48
David Mayerich
bug fixes for Net...
|
96
97
|
///adds a magnitude to each point in the cylinder
|
7115e973
David Mayerich
removed the magni...
|
98
|
/*void add_mag(V val = 0) {
|
eb554b48
David Mayerich
bug fixes for Net...
|
99
100
|
if (M.size() == 0) M.resize(size()); //if the magnitude vector isn't initialized, resize it to match the centerline
for (size_t i = 0; i < size(); i++) //for each point
|
7115e973
David Mayerich
removed the magni...
|
101
102
|
R[i].push_back(val); //add this value to the magnitude vector at each point
}*/
|
eb554b48
David Mayerich
bug fixes for Net...
|
103
104
|
//adds a magnitude based on a list of magnitudes for each point
|
7115e973
David Mayerich
removed the magni...
|
105
|
/*void add_mag(std::vector<T> val) {
|
eb554b48
David Mayerich
bug fixes for Net...
|
106
107
|
if (M.size() == 0) M.resize(size()); //if the magnitude vector isn't initialized, resize it to match the centerline
for (size_t i = 0; i < size(); i++) //for each point
|
7115e973
David Mayerich
removed the magni...
|
108
109
|
R[i].push_back(val[i]); //add this value to the magnitude vector at each point
}*/
|
eb554b48
David Mayerich
bug fixes for Net...
|
110
111
|
//sets the value of magnitude m at point i
|
7115e973
David Mayerich
removed the magni...
|
112
113
|
void set_r(size_t i, T r) {
R[i] = r;
|
eb554b48
David Mayerich
bug fixes for Net...
|
114
115
|
}
|
7115e973
David Mayerich
removed the magni...
|
116
|
/*size_t nmags() {
|
eb554b48
David Mayerich
bug fixes for Net...
|
117
|
if (M.size() == 0) return 0;
|
7115e973
David Mayerich
removed the magni...
|
118
119
|
else return R[0].size();
}*/
|
eb554b48
David Mayerich
bug fixes for Net...
|
120
|
|
e21d1051
David Mayerich
cylinder can buil...
|
121
|
///Returns a circle representing the cylinder cross section at point i
|
7115e973
David Mayerich
removed the magni...
|
122
123
|
stim::circle<T> circ(size_t i) {
return stim::circle<T>(at(i), R[i], d(i), U[i]);
|
e21d1051
David Mayerich
cylinder can buil...
|
124
125
126
|
}
///Returns an OBJ object representing the cylinder with a radial tesselation value of N using magnitude m
|
7115e973
David Mayerich
removed the magni...
|
127
|
stim::obj<T> OBJ(size_t N) {
|
e21d1051
David Mayerich
cylinder can buil...
|
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
|
stim::obj<T> out; //create an OBJ object
stim::circle<T> c0, c1;
std::vector< stim::vec3<T> > p0, p1; //allocate space for the point sets representing the circles bounding each cylinder segment
T u0, u1, v0, v1; //allocate variables to store running texture coordinates
for (size_t i = 1; i < size(); i++) { //for each line segment in the cylinder
c0 = circ(i - 1); //get the two circles bounding the line segment
c1 = circ(i);
p0 = c0.points(N); //get t points for each of the end caps
p1 = c1.points(N);
u0 = L[i - 1] / length(); //calculate the texture coordinate (u, v) where u runs along the cylinder length
u1 = L[i] / length();
for (size_t n = 1; n < N; n++) { //for each point in the circle
v0 = (double)(n-1) / (double)(N - 1); //v texture coordinate runs around the cylinder
v1 = (double)(n) / (double)(N - 1);
out.Begin(OBJ_FACE); //start a face (quad)
out.TexCoord(u0, v0);
out.Vertex(p0[n - 1][0], p0[n - 1][1], p0[n - 1][2]); //output the points composing a strip of quads wrapping the cylinder segment
out.TexCoord(u1, v0);
out.Vertex(p1[n - 1][0], p1[n - 1][1], p1[n - 1][2]);
out.TexCoord(u0, v1);
out.Vertex(p1[n][0], p1[n][1], p1[n][2]);
out.TexCoord(u1, v1);
out.Vertex(p0[n][0], p0[n][1], p0[n][2]);
out.End();
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
156
|
}
|
e21d1051
David Mayerich
cylinder can buil...
|
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
|
v0 = (double)(N - 2) / (double)(N - 1); //v texture coordinate runs around the cylinder
v1 = 1.0;
out.Begin(OBJ_FACE);
out.TexCoord(u0, v0);
out.Vertex(p0[N - 1][0], p0[N - 1][1], p0[N - 1][2]); //output the points composing a strip of quads wrapping the cylinder segment
out.TexCoord(u1, v0);
out.Vertex(p1[N - 1][0], p1[N - 1][1], p1[N - 1][2]);
out.TexCoord(u0, v1);
out.Vertex(p1[0][0], p1[0][1], p1[0][2]);
out.TexCoord(u1, v1);
out.Vertex(p0[0][0], p0[0][1], p0[0][2]);
out.End();
}
return out;
}
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
173
|
|
eb554b48
David Mayerich
bug fixes for Net...
|
174
175
176
177
178
|
std::string str() {
std::stringstream ss;
size_t N = std::vector< stim::vec3<T> >::size();
ss << "---------[" << N << "]---------" << std::endl;
for (size_t i = 0; i < N; i++)
|
7115e973
David Mayerich
removed the magni...
|
179
|
ss << std::vector< stim::vec3<T> >::at(i) << " r = " << R[i] << " u = " << U[i] << std::endl;
|
eb554b48
David Mayerich
bug fixes for Net...
|
180
181
182
183
184
185
186
|
ss << "--------------------" << std::endl;
return ss.str();
}
/// Integrates a magnitude value along the cylinder.
/// @param m is the magnitude value to be integrated (this is usually the radius)
|
7115e973
David Mayerich
removed the magni...
|
187
|
T integrate() {
|
58674036
Jiaming Guo
fix minor error f...
|
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
|
T sum = 0; //initialize the integral to zero
if (L.size() == 1)
return sum;
else {
T m0, m1; //allocate space for both magnitudes in a single segment
m0 = R[0]; //initialize the first point and magnitude to the first point in the cylinder
T len = L[1]; //allocate space for the segment length
for (unsigned p = 1; p < size(); p++) { //for every consecutive point in the cylinder
m1 = R[p];
if (p > 1) len = (L[p] - L[p - 1]); //calculate the segment length using the L array
sum += (m0 + m1) / (T)2.0 * len; //add the average magnitude, weighted by the segment length
m0 = m1; //move to the next segment by shifting points
}
return sum; //return the integral
|
eb554b48
David Mayerich
bug fixes for Net...
|
204
|
}
|
eb554b48
David Mayerich
bug fixes for Net...
|
205
206
207
208
209
210
211
212
|
}
/// 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) {
cylinder<T> c = stim::centerline<T>::resample(spacing); //resample the centerline and use it to create a new cylinder
|
7115e973
David Mayerich
removed the magni...
|
213
214
215
216
217
|
//size_t nm = nmags(); //get the number of magnitude values in the current cylinder
//if (nm > 0) { //if there are magnitude values
// std::vector<T> magvec(nm, 0); //create a magnitude vector for a single point
// c.M.resize(c.size(), magvec); //allocate space for a magnitude vector at each point of the new cylinder
//}
|
eb554b48
David Mayerich
bug fixes for Net...
|
218
219
220
221
222
|
T l, t;
for (size_t i = 0; i < c.size(); i++) { //for each point in the new cylinder
l = c.L[i]; //get the length along the new cylinder
t = l / length(); //calculate the parameter value along the new cylinder
|
7115e973
David Mayerich
removed the magni...
|
223
224
225
|
//for (size_t mag = 0; mag < nm; mag++) { //for each magnitude value
c.R[i] = r(t); //retrieve the interpolated magnitude from the current cylinder and store it in the new one
//}
|
eb554b48
David Mayerich
bug fixes for Net...
|
226
|
}
|
0d0ef1a9
David Mayerich
Adapted classes t...
|
227
228
|
return c;
}
|
eb554b48
David Mayerich
bug fixes for Net...
|
229
|
|
0d0ef1a9
David Mayerich
Adapted classes t...
|
230
231
232
233
234
235
236
237
238
239
|
std::vector< stim::cylinder<T> > split(unsigned int idx) {
unsigned N = size();
std::vector< stim::centerline<T> > LL;
LL.resize(2);
LL = (*this).centerline<T>::split(idx);
std::vector< stim::cylinder<T> > C(LL.size());
unsigned i = 0;
C[0] = LL[0];
//C[0].R.resize(idx);
|
da89bce9
Jiaming Guo
fixed splitting bugs
|
240
|
for (; i < idx + 1; i++) {
|
0d0ef1a9
David Mayerich
Adapted classes t...
|
241
242
243
244
245
246
247
|
//for(unsigned d = 0; d < 3; d++)
//C[0][i][d] = LL[0].c[i][d];
C[0].R[i] = R[i];
//C[0].R[i].resize(1);
}
if (C.size() == 2) {
C[1] = LL[1];
|
da89bce9
Jiaming Guo
fixed splitting bugs
|
248
|
i--;
|
0d0ef1a9
David Mayerich
Adapted classes t...
|
249
250
251
252
253
254
255
256
|
//C[1].M.resize(N - idx);
for (; i < N; i++) {
//for(unsigned d = 0; d < 3; d++)
//C[1][i - idx][d] = LL[1].c[i - idx][d];
C[1].R[i - idx] = R[i];
//C[1].M[i - idx].resize(1);
}
}
|
eb554b48
David Mayerich
bug fixes for Net...
|
257
|
|
0d0ef1a9
David Mayerich
Adapted classes t...
|
258
|
return C;
|
eb554b48
David Mayerich
bug fixes for Net...
|
259
260
|
}
|
e21d1051
David Mayerich
cylinder can buil...
|
261
262
263
|
/*
///inits the cylinder from a list of points (std::vector of stim::vec3 --inP) and magnitudes (inM)
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
264
|
void
|
e21d1051
David Mayerich
cylinder can buil...
|
265
266
267
|
init(centerline inP, std::vector< std::vector<T> > inM) {
M = inM; //the magnitude vector can be copied directly
(*this) = inP; //the centerline can be copied to this class directly
|
308a743c
David Mayerich
fixed class compa...
|
268
269
|
stim::vec3<float> v1;
stim::vec3<float> v2;
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
270
|
e.resize(inP.size());
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
271
272
273
274
|
norms.resize(inP.size());
Us.resize(inP.size());
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
275
276
|
if(inP.size() < 2)
return;
|
8495a970
Pavel Govyadinov
added comments to...
|
277
278
|
//calculate each L.
|
1b2858a2
David Mayerich
fixed some commen...
|
279
280
281
|
L.resize(inP.size()); //the number of precomputed lengths will equal the number of points
T temp = (T)0; //length up to that point
L[0] = temp;
|
308a743c
David Mayerich
fixed class compa...
|
282
|
for(size_t i = 1; i < L.size(); i++)
|
10654a1f
Pavel Govyadinov
added the necessa...
|
283
|
{
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
284
|
temp += (inP[i-1] - inP[i]).len();
|
bf23ee36
Pavel Govyadinov
more minor bug fi...
|
285
|
L[i] = temp;
|
10654a1f
Pavel Govyadinov
added the necessa...
|
286
|
}
|
9c97e126
David Mayerich
added an axis-ali...
|
287
|
|
308a743c
David Mayerich
fixed class compa...
|
288
|
stim::vec3<T> dr = (inP[1] - inP[0]).norm();
|
7115e973
David Mayerich
removed the magni...
|
289
|
s = stim::circle<T>(inP[0], inR[0][0], dr, stim::vec3<T>(1,0,0));
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
290
|
norms[0] = s.N;
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
291
|
e[0] = s;
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
292
|
Us[0] = e[0].U;
|
308a743c
David Mayerich
fixed class compa...
|
293
|
for(size_t i = 1; i < inP.size()-1; i++)
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
294
295
296
297
298
299
|
{
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);
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
300
301
302
|
norms[i] = s.N;
|
7115e973
David Mayerich
removed the magni...
|
303
|
s.scale(inR[i][0]/inR[i-1][0]);
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
304
|
e[i] = s;
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
305
|
Us[i] = e[i].U;
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
306
307
308
309
310
311
|
}
int j = inP.size()-1;
s.center(inP[j]);
dr = (inP[j] - inP[j-1]).norm();
s.normal(dr);
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
312
313
314
|
norms[j] = s.N;
|
7115e973
David Mayerich
removed the magni...
|
315
|
s.scale(inR[j][0]/inR[j-1][0]);
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
316
|
e[j] = s;
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
317
|
Us[j] = e[j].U;
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
318
319
|
}
|
8495a970
Pavel Govyadinov
added comments to...
|
320
|
///returns the direction vector at point idx.
|
308a743c
David Mayerich
fixed class compa...
|
321
|
stim::vec3<T>
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
322
323
324
325
|
d(int idx)
{
if(idx == 0)
{
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
326
327
328
329
330
331
332
|
stim::vec3<T> temp(
c[idx+1][0]-c[idx][0],
c[idx+1][1]-c[idx][1],
c[idx+1][2]-c[idx][2]
);
// return (e[idx+1].P - e[idx].P).norm();
return (temp.norm());
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
333
|
}
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
334
|
else if(idx == N-1)
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
335
|
{
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
336
337
338
339
340
341
342
|
stim::vec3<T> temp(
c[idx][0]-c[idx+1][0],
c[idx][1]-c[idx+1][1],
c[idx][2]-c[idx+1][2]
);
// return (e[idx].P - e[idx-1].P).norm();
return (temp.norm());
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
343
344
345
|
}
else
{
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
346
|
// return (e[idx+1].P - e[idx].P).norm();
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
|
// stim::vec3<float> v1 = (e[idx].P-e[idx-1].P).norm();
stim::vec3<T> v1(
c[idx][0]-c[idx-1][0],
c[idx][1]-c[idx-1][1],
c[idx][2]-c[idx-1][2]
);
// stim::vec3<float> v2 = (e[idx+1].P-e[idx].P).norm();
stim::vec3<T> v2(
c[idx+1][0]-c[idx][0],
c[idx+1][1]-c[idx][1],
c[idx+1][2]-c[idx][2]
);
return (v1.norm()+v2.norm()).norm();
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
362
363
364
|
}
// return e[idx].N;
|
10654a1f
Pavel Govyadinov
added the necessa...
|
365
366
|
}
|
308a743c
David Mayerich
fixed class compa...
|
367
|
stim::vec3<T>
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
368
|
d(T l, int idx)
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
369
|
{
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
370
|
if(idx == 0 || idx == N-1)
|
10654a1f
Pavel Govyadinov
added the necessa...
|
371
|
{
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
372
373
|
return norms[idx];
// return e[idx].N;
|
10654a1f
Pavel Govyadinov
added the necessa...
|
374
|
}
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
375
376
|
else
{
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
377
|
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
378
|
T rat = (l-L[idx])/(L[idx+1]-L[idx]);
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
379
380
|
return( norms[idx] + (norms[idx+1] - norms[idx])*rat);
// return( e[idx].N + (e[idx+1].N - e[idx].N)*rat);
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
381
|
}
|
10654a1f
Pavel Govyadinov
added the necessa...
|
382
383
|
}
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
384
|
|
8495a970
Pavel Govyadinov
added comments to...
|
385
386
|
///finds the index of the point closest to the length l on the lower bound.
///binary search.
|
10654a1f
Pavel Govyadinov
added the necessa...
|
387
|
int
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
388
389
390
391
392
|
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...
|
393
|
while(i > 0 && i < L.size()-1)
|
10654a1f
Pavel Govyadinov
added the necessa...
|
394
|
{
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
395
396
|
// std::cerr << "Trying " << i << std::endl;
// std::cerr << "l is " << l << ", L[" << i << "]" << L[i] << std::endl;
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
397
|
if(l < L[i])
|
10654a1f
Pavel Govyadinov
added the necessa...
|
398
|
{
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
399
|
max = i;
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
400
|
i = min+(max-min)/2;
|
10654a1f
Pavel Govyadinov
added the necessa...
|
401
|
}
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
402
|
else if(L[i] <= l && L[i+1] >= l)
|
10654a1f
Pavel Govyadinov
added the necessa...
|
403
404
405
406
407
|
{
break;
}
else
{
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
408
|
min = i;
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
409
|
i = min+(max-min)/2;
|
10654a1f
Pavel Govyadinov
added the necessa...
|
410
411
412
413
|
}
}
return i;
}
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
414
|
|
9c97e126
David Mayerich
added an axis-ali...
|
415
416
|
public:
///default constructor
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
417
|
cylinder()
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
418
|
// : centerline<T>()
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
419
420
421
|
{
}
|
9c97e126
David Mayerich
added an axis-ali...
|
422
|
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
423
|
///constructor to create a cylinder from a set of points, radii, and the number of sides for the cylinder.
|
7d3162a2
Pavel Govyadinov
fixed majority of...
|
424
|
///@param inP: Vector of stim vec3 composing the points of the centerline.
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
425
|
///@param inM: Vector of stim vecs composing the radii of the centerline.
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
426
427
428
|
cylinder(std::vector<stim::vec3<T> > inP, std::vector<stim::vec<T> > inM)
: centerline<T>(inP)
{
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
429
430
431
|
init(inP, inM);
}
|
09049866
Pavel Govyadinov
fixed a drawing b...
|
432
433
434
435
436
437
438
|
///constructor to create a cylinder from a set of points, radii, and the number of sides for the cylinder.
///@param inP: Vector of stim vec3 composing the points of the centerline.
///@param inM: Vector of stim vecs composing the radii of the centerline.
cylinder(std::vector<stim::vec3<T> > inP, std::vector< T > inM)
: centerline<T>(inP)
{
std::vector<stim::vec<T> > temp;
|
a36c3850
Pavel Govyadinov
fixed the bug tha...
|
439
|
stim::vec<T> zero(0.0,0.0);
|
09049866
Pavel Govyadinov
fixed a drawing b...
|
440
|
temp.resize(inM.size(), zero);
|
47fe6300
Pavel Govyadinov
fixed an accident...
|
441
|
for(int i = 0; i < inM.size(); i++)
|
7115e973
David Mayerich
removed the magni...
|
442
|
temp[i][0] = inR[i];
|
09049866
Pavel Govyadinov
fixed a drawing b...
|
443
444
445
|
init(inP, temp);
}
|
7d3162a2
Pavel Govyadinov
fixed majority of...
|
446
|
|
58ee2b21
David Mayerich
incorporated stim...
|
447
|
///Constructor defines a cylinder with centerline inP and magnitudes of zero
|
7d3162a2
Pavel Govyadinov
fixed majority of...
|
448
|
///@param inP: Vector of stim vec3 composing the points of the centerline
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
449
450
451
|
cylinder(std::vector< stim::vec3<T> > inP)
: centerline<T>(inP)
{
|
96e8ab1f
David Mayerich
fixed stim::cylin...
|
452
|
std::vector< stim::vec<T> > inM; //create an array of arbitrary magnitudes
|
9c97e126
David Mayerich
added an axis-ali...
|
453
454
455
456
457
|
stim::vec<T> zero;
zero.push_back(0);
inM.resize(inP.size(), zero); //initialize the magnitude values to zero
|
58ee2b21
David Mayerich
incorporated stim...
|
458
459
460
|
init(inP, inM);
}
|
e21d1051
David Mayerich
cylinder can buil...
|
461
462
463
464
465
|
//assignment operator creates a cylinder from a centerline (default radius is zero)
cylinder& operator=(stim::centerline<T> c) {
init(c);
}
|
58ee2b21
David Mayerich
incorporated stim...
|
466
467
468
|
///Returns the number of points on the cylinder centerline
unsigned int size(){
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
469
|
return N;
|
58ee2b21
David Mayerich
incorporated stim...
|
470
471
|
}
|
e21d1051
David Mayerich
cylinder can buil...
|
472
|
|
10654a1f
Pavel Govyadinov
added the necessa...
|
473
|
///Returns a position vector at the given p-value (p value ranges from 0 to 1).
|
8495a970
Pavel Govyadinov
added comments to...
|
474
475
|
///interpolates the position along the line.
///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
|
308a743c
David Mayerich
fixed class compa...
|
476
|
stim::vec3<T>
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
477
478
|
p(T pvalue)
{
|
bf23ee36
Pavel Govyadinov
more minor bug fi...
|
479
|
if(pvalue < 0.0 || pvalue > 1.0)
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
480
|
{
|
308a743c
David Mayerich
fixed class compa...
|
481
|
return stim::vec3<float>(-1,-1,-1);
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
482
|
}
|
10654a1f
Pavel Govyadinov
added the necessa...
|
483
484
|
T l = pvalue*L[L.size()-1];
int idx = findIdx(l);
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
485
|
return (p(l,idx));
|
e21d1051
David Mayerich
cylinder can buil...
|
486
|
}
|
10654a1f
Pavel Govyadinov
added the necessa...
|
487
|
|
8495a970
Pavel Govyadinov
added comments to...
|
488
489
490
491
|
///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.
|
308a743c
David Mayerich
fixed class compa...
|
492
|
stim::vec3<T>
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
493
494
495
|
p(T l, int idx)
{
T rat = (l-L[idx])/(L[idx+1]-L[idx]);
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
496
497
498
499
500
501
502
503
504
505
506
507
508
|
stim::vec3<T> v1(
c[idx][0],
c[idx][1],
c[idx][2]
);
stim::vec3<T> v2(
c[idx+1][0],
c[idx+1][1],
c[idx+1][2]
);
// return( e[idx].P + (e[idx+1].P-e[idx].P)*rat);
return( v1 + (v2-v1)*rat);
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
509
510
|
// return(
// return (pos[idx] + (pos[idx+1]-pos[idx])*((l-L[idx])/(L[idx+1]- L[idx])));
|
10654a1f
Pavel Govyadinov
added the necessa...
|
511
512
513
|
}
///Returns a radius at the given p-value (p value ranges from 0 to 1).
|
8495a970
Pavel Govyadinov
added comments to...
|
514
515
|
///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...
|
516
|
T
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
517
518
|
r(T pvalue)
{
|
98eecaa9
David Mayerich
VS and win32 updates
|
519
520
521
522
|
if(pvalue < 0.0 || pvalue > 1.0){
std::cerr<<"Error, value "<<pvalue<<" is outside of [0 1]."<<std::endl;
exit(1);
}
|
10654a1f
Pavel Govyadinov
added the necessa...
|
523
524
|
T l = pvalue*L[L.size()-1];
int idx = findIdx(l);
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
525
|
return (r(l,idx));
|
e21d1051
David Mayerich
cylinder can buil...
|
526
|
}
|
10654a1f
Pavel Govyadinov
added the necessa...
|
527
|
|
8495a970
Pavel Govyadinov
added comments to...
|
528
529
530
531
|
///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...
|
532
|
T
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
533
534
535
|
r(T l, int idx)
{
T rat = (l-L[idx])/(L[idx+1]-L[idx]);
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
536
537
538
539
540
541
542
543
|
T v1 = (e[idx].U.len() + (e[idx+1].U.len() - e[idx].U.len())*rat);
T v3 = (Us[idx].len() + (Us[idx+1].len() - Us[idx].len())*rat);
T v2 = (mags[idx][0] + (mags[idx+1][0]-mags[idx][0])*rat);
// std::cout << (float)v1 = (float) v2 << std::endl;
if(abs(v3 - v1) >= 10e-6)
{
std::cout << "-------------------------" << std::endl;
std::cout << e[idx].str() << std::endl << std::endl;
|
b50c840e
David Mayerich
eliminated ANN fr...
|
544
|
std::cout << Us[idx].str() << std::endl;
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
545
546
547
548
549
|
std::cout << (float)v1 - (float) v2 << std::endl;
std::cout << "failed" << std::endl;
}
// std::cout << e[idx].U.len() << " " << mags[idx][0] << std::endl;
// std::cout << v2 << std::endl;
|
24dee661
Pavel Govyadinov
finished the debu...
|
550
|
return(v2);
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
551
552
|
// return (mags[idx][0] + (mags[idx+1][0]-mags[idx][0])*rat);
// (
|
66fe63f3
Pavel Govyadinov
fixed minor typo
|
553
|
}
|
10654a1f
Pavel Govyadinov
added the necessa...
|
554
|
|
9c97e126
David Mayerich
added an axis-ali...
|
555
556
557
558
559
560
561
562
563
564
|
/// 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){
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
565
|
for(unsigned int p = 0; p < N; p++)
|
9c97e126
David Mayerich
added an axis-ali...
|
566
567
568
569
570
571
572
573
574
575
576
|
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...
|
577
578
579
580
581
|
/// Returns the number of magnitude values at each point
unsigned nmags(){
return mags[0].size();
}
|
10654a1f
Pavel Govyadinov
added the necessa...
|
582
|
///returns the position of the point with a given pvalue and theta on the surface
|
8495a970
Pavel Govyadinov
added comments to...
|
583
584
585
|
///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.
|
308a743c
David Mayerich
fixed class compa...
|
586
|
stim::vec3<T>
|
10654a1f
Pavel Govyadinov
added the necessa...
|
587
588
|
surf(T pvalue, T theta)
{
|
bf23ee36
Pavel Govyadinov
more minor bug fi...
|
589
|
if(pvalue < 0.0 || pvalue > 1.0)
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
590
|
{
|
308a743c
David Mayerich
fixed class compa...
|
591
|
return stim::vec3<float>(-1,-1,-1);
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
592
|
} else {
|
10654a1f
Pavel Govyadinov
added the necessa...
|
593
594
|
T l = pvalue*L[L.size()-1];
int idx = findIdx(l);
|
308a743c
David Mayerich
fixed class compa...
|
595
|
stim::vec3<T> ps = p(l, idx);
|
10654a1f
Pavel Govyadinov
added the necessa...
|
596
|
T m = r(l, idx);
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
597
598
599
600
|
s = e[idx];
s.center(ps);
s.normal(d(l, idx));
s.scale(m/e[idx].U.len());
|
10654a1f
Pavel Govyadinov
added the necessa...
|
601
|
return(s.p(theta));
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
602
|
}
|
10654a1f
Pavel Govyadinov
added the necessa...
|
603
604
|
}
|
8495a970
Pavel Govyadinov
added comments to...
|
605
|
///returns a vector of points necessary to create a circle at every position in the fiber.
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
606
|
///@param sides: the number of sides of each circle.
|
308a743c
David Mayerich
fixed class compa...
|
607
|
std::vector<std::vector<vec3<T> > >
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
608
609
|
getPoints(int sides)
{
|
308a743c
David Mayerich
fixed class compa...
|
610
|
std::vector<std::vector <vec3<T> > > points;
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
611
612
|
points.resize(N);
for(int i = 0; i < N; i++)
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
613
|
{
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
614
|
points[i] = e[i].getPoints(sides);
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
615
|
}
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
616
|
return points;
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
617
|
}
|
58ee2b21
David Mayerich
incorporated stim...
|
618
|
|
8c4f5d84
Pavel Govyadinov
fixed the issue w...
|
619
620
621
622
623
624
|
///returns the total length of the line at index j.
T
getl(int j)
{
return (L[j]);
}
|
58ee2b21
David Mayerich
incorporated stim...
|
625
626
|
/// Allows a point on the centerline to be accessed using bracket notation
|
308a743c
David Mayerich
fixed class compa...
|
627
|
vec3<T> operator[](unsigned int i){
|
3e5d3ad3
Pavel Govyadinov
merged the change...
|
628
|
return e[i].P;
|
58ee2b21
David Mayerich
incorporated stim...
|
629
630
631
632
633
634
635
|
}
/// Returns the total length of the cylinder centerline
T length(){
return L.back();
}
|
9c97e126
David Mayerich
added an axis-ali...
|
636
637
638
639
640
641
642
|
/// 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
|
308a743c
David Mayerich
fixed class compa...
|
643
|
//vec3<T> p0, p1; //allocate space for both points in a single segment
|
9c97e126
David Mayerich
added an axis-ali...
|
644
645
646
647
648
649
650
|
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
|
0311ab81
Pavel Govyadinov
fixed some issues...
|
651
|
for(unsigned p = 1; p < N; p++){
|
9c97e126
David Mayerich
added an axis-ali...
|
652
653
654
655
656
657
658
|
//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
|
308a743c
David Mayerich
fixed class compa...
|
659
|
M += (m0 + m1)/(T)2.0 * len;
|
9c97e126
David Mayerich
added an axis-ali...
|
660
661
662
663
664
665
666
667
668
669
670
671
672
673
|
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...
|
674
675
676
677
678
|
/// 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){
|
308a743c
David Mayerich
fixed class compa...
|
679
|
std::vector< vec3<T> > result;
|
9c97e126
David Mayerich
added an axis-ali...
|
680
|
|
1b2858a2
David Mayerich
fixed some commen...
|
681
|
vec3<T> p0 = e[0].P; //initialize p0 to the first point on the centerline
|
308a743c
David Mayerich
fixed class compa...
|
682
|
vec3<T> p1;
|
1b2858a2
David Mayerich
fixed some commen...
|
683
|
unsigned N = size(); //number of points in the current centerline
|
58ee2b21
David Mayerich
incorporated stim...
|
684
685
686
|
//for each line segment on the centerline
for(unsigned int i = 1; i < N; i++){
|
1b2858a2
David Mayerich
fixed some commen...
|
687
|
p1 = e[i].P; //get the second point in the line segment
|
58ee2b21
David Mayerich
incorporated stim...
|
688
|
|
1b2858a2
David Mayerich
fixed some commen...
|
689
690
|
vec3<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)
|
58ee2b21
David Mayerich
incorporated stim...
|
691
|
|
308a743c
David Mayerich
fixed class compa...
|
692
|
size_t nsteps = (size_t)std::ceil(d / spacing); //calculate the number of steps to take along the segment to meet the spacing criteria
|
1b2858a2
David Mayerich
fixed some commen...
|
693
|
T stepsize = (T)1.0 / nsteps; //calculate the parametric step size between new centerline points
|
58ee2b21
David Mayerich
incorporated stim...
|
694
695
696
|
//for each step along the line segment
for(unsigned s = 0; s < nsteps; s++){
|
1b2858a2
David Mayerich
fixed some commen...
|
697
698
|
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
|
58ee2b21
David Mayerich
incorporated stim...
|
699
700
|
}
|
1b2858a2
David Mayerich
fixed some commen...
|
701
|
p0 = p1; //shift the points to move to the next line segment
|
58ee2b21
David Mayerich
incorporated stim...
|
702
703
|
}
|
1b2858a2
David Mayerich
fixed some commen...
|
704
|
result.push_back(e[size() - 1].P); //push the last point in the centerline
|
9c97e126
David Mayerich
added an axis-ali...
|
705
|
|
58ee2b21
David Mayerich
incorporated stim...
|
706
707
|
return cylinder<T>(result);
|
e21d1051
David Mayerich
cylinder can buil...
|
708
|
}*/
|
9c97e126
David Mayerich
added an axis-ali...
|
709
|
|
26aee9ed
Pavel Govyadinov
changed circle cl...
|
710
|
|
5b6c317a
Pavel Govyadinov
implemented the c...
|
711
712
713
|
};
}
|
c72184d1
Jiaming Guo
add many function...
|
714
|
#endif
|