Blame view

stim/visualization/cylinder.h 10.2 KB
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
  			init(inP, inM);
  		}
  
  
  		///Returns the number of points on the cylinder centerline
  
  		unsigned int size(){
9b766f1f   Pavel Govyadinov   completed merge f...
167
  			return e.size();
58ee2b21   David Mayerich   incorporated stim...
168
169
  		}
  
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
  		///	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){
9b766f1f   Pavel Govyadinov   completed merge f...
234
  			for(unsigned int p = 0; p < e.size(); p++)
9c97e126   David Mayerich   added an axis-ali...
235
236
237
238
239
240
241
242
243
244
245
  				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
  		/// 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
9b766f1f   Pavel Govyadinov   completed merge f...
320
  			for(unsigned p = 1; p < e.size(); p++){
9c97e126   David Mayerich   added an axis-ali...
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
  
  				//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
  
9b766f1f   Pavel Govyadinov   completed merge f...
350
  			vec<T> p0 = e[0].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