Branch data Line data Source code
1 : : /*=========================================================================
2 : :
3 : : Module: $RCSfile: VerdictVector.hpp,v $
4 : :
5 : : Copyright (c) 2006 Sandia Corporation.
6 : : All rights reserved.
7 : : See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
8 : :
9 : : This software is distributed WITHOUT ANY WARRANTY; without even
10 : : the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 : : PURPOSE. See the above copyright notice for more information.
12 : :
13 : : =========================================================================*/
14 : :
15 : : /*
16 : : *
17 : : * VerdictVector.hpp contains declarations of vector operations
18 : : *
19 : : * This file is part of VERDICT
20 : : *
21 : : */
22 : :
23 : : #ifndef VERDICTVECTOR_HPP
24 : : #define VERDICTVECTOR_HPP
25 : :
26 : : #include "moab/verdict.h"
27 : : #include <assert.h>
28 : : #include <math.h>
29 : :
30 : : class VerdictVector;
31 : : typedef void ( VerdictVector::*transform_function )( double gamma, double gamma2 );
32 : : // a pointer to some function that transforms the point,
33 : : // taking a double parameter. e.g. blow_out, rotate, and scale_angle
34 : :
35 : : class VerdictVector
36 : : {
37 : : public:
38 : : //- Heading: Constructors and Destructor
39 : : VerdictVector(); //- Default constructor.
40 : :
41 : : VerdictVector( const double x, const double y, const double z );
42 : : //- Constructor: create vector from three components
43 : :
44 : : VerdictVector( const double xyz[3] );
45 : : //- Constructor: create vector from tuple
46 : :
47 : : VerdictVector( const VerdictVector& tail, const VerdictVector& head );
48 : : //- Constructor for a VerdictVector starting at tail and pointing
49 : : //- to head.
50 : :
51 : : VerdictVector( const VerdictVector& copy_from ); //- Copy Constructor
52 : :
53 : : //- Heading: Set and Inquire Functions
54 : : void set( const double xv, const double yv, const double zv );
55 : : //- Change vector components to {x}, {y}, and {z}
56 : :
57 : : void set( const double xyz[3] );
58 : : //- Change vector components to xyz[0], xyz[1], xyz[2]
59 : :
60 : : void set( const VerdictVector& tail, const VerdictVector& head );
61 : : //- Change vector to go from tail to head.
62 : :
63 : : void set( const VerdictVector& to_copy );
64 : : //- Same as operator=(const VerdictVector&)
65 : :
66 : : double x() const; //- Return x component of vector
67 : :
68 : : double y() const; //- Return y component of vector
69 : :
70 : : double z() const; //- Return z component of vector
71 : :
72 : : void get_xyz( double& x, double& y, double& z ); //- Get x, y, z components
73 : : void get_xyz( double xyz[3] ); //- Get xyz tuple
74 : :
75 : : double& r(); //- Return r component of vector, if (r,theta) format
76 : :
77 : : double& theta(); //- Return theta component of vector, if (r,theta) format
78 : :
79 : : void x( const double xv ); //- Set x component of vector
80 : :
81 : : void y( const double yv ); //- Set y component of vector
82 : :
83 : : void z( const double zv ); //- Set z component of vector
84 : :
85 : : void r( const double xv ); //- Set r component of vector, if (r,theta) format
86 : :
87 : : void theta( const double yv ); //- Set theta component of vector, if (r,theta) format
88 : :
89 : : void xy_to_rtheta();
90 : : //- convert from cartesian to polar coordinates, just 2d for now
91 : : //- theta is in [0,2 PI)
92 : :
93 : : void rtheta_to_xy();
94 : : //- convert from polar to cartesian coordinates, just 2d for now
95 : :
96 : : void scale_angle( double gamma, double );
97 : : //- tranform_function.
98 : : //- transform (x,y) to (r,theta) to (r,gamma*theta) to (x',y')
99 : : //- plus some additional scaling so long chords won't cross short ones
100 : :
101 : : void blow_out( double gamma, double gamma2 = 0.0 );
102 : : //- transform_function
103 : : //- blow points radially away from the origin,
104 : :
105 : : void rotate( double angle, double );
106 : : //- transform function.
107 : : //- transform (x,y) to (r,theta) to (r,theta+angle) to (x',y')
108 : :
109 : : void reflect_about_xaxis( double dummy, double );
110 : : //- dummy argument to make it a transform function
111 : :
112 : : double normalize();
113 : : //- Normalize (set magnitude equal to 1) vector - return the magnitude
114 : :
115 : : VerdictVector& length( const double new_length );
116 : : //- Change length of vector to {new_length}. Can be used to move a
117 : : //- location a specified distance from the origin in the current
118 : : //- orientation.
119 : :
120 : : double length() const;
121 : : //- Calculate the length of the vector.
122 : : //- Use {length_squared()} if only comparing lengths, not adding.
123 : :
124 : : double distance_between( const VerdictVector& test_vector );
125 : : //- Calculate the distance from the head of one vector
126 : : // to the head of the test_vector.
127 : :
128 : : double length_squared() const;
129 : : //- Calculate the squared length of the vector.
130 : : //- Faster than {length()} since it eliminates the square root if
131 : : //- only comparing other lengths.
132 : :
133 : : double interior_angle( const VerdictVector& otherVector );
134 : : //- Calculate the interior angle: acos((a%b)/(|a||b|))
135 : : //- Returns angle in degrees.
136 : :
137 : : double vector_angle_quick( const VerdictVector& vec1, const VerdictVector& vec2 );
138 : : //- Calculate the interior angle between the projections of
139 : : //- {vec1} and {vec2} onto the plane defined by the {this} vector.
140 : : //- The angle returned is the right-handed angle around the {this}
141 : : //- vector from {vec1} to {vec2}. Angle is in radians.
142 : :
143 : : double vector_angle( const VerdictVector& vector1, const VerdictVector& vector2 ) const;
144 : : //- Compute the angle between the projections of {vector1} and {vector2}
145 : : //- onto the plane defined by *this. The angle is the
146 : : //- right-hand angle, in radians, about *this from {vector1} to {vector2}.
147 : :
148 : : void perpendicular_z();
149 : : //- Transform this vector to a perpendicular one, leaving
150 : : //- z-component alone. Rotates clockwise about the z-axis by pi/2.
151 : :
152 : : void print_me();
153 : : //- Prints out the coordinates of this vector.
154 : :
155 : : void orthogonal_vectors( VerdictVector& vector2, VerdictVector& vector3 );
156 : : //- Finds 2 (arbitrary) vectors that are orthogonal to this one
157 : :
158 : : void next_point( const VerdictVector& direction, double distance, VerdictVector& out_point );
159 : : //- Finds the next point in space based on *this* point (starting point),
160 : : //- a direction and the distance to extend in the direction. The
161 : : //- direction vector need not be a unit vector. The out_point can be
162 : : //- "*this" (i.e., modify point in place).
163 : :
164 : : bool within_tolerance( const VerdictVector& vectorPtr2, double tolerance ) const;
165 : : //- Compare two vectors to see if they are spatially equal. They
166 : : //- compare if x, y, and z are all within tolerance.
167 : :
168 : : //- Heading: Operator Overloads *****************************
169 : : VerdictVector& operator+=( const VerdictVector& vec );
170 : : //- Compound Assignment: addition: {this = this + vec}
171 : :
172 : : VerdictVector& operator-=( const VerdictVector& vec );
173 : : //- Compound Assignment: subtraction: {this = this - vec}
174 : :
175 : : VerdictVector& operator*=( const VerdictVector& vec );
176 : : //- Compound Assignment: cross product: {this = this * vec},
177 : : //- non-commutative
178 : :
179 : : VerdictVector& operator*=( const double scalar );
180 : : //- Compound Assignment: multiplication: {this = this * scalar}
181 : :
182 : : VerdictVector& operator/=( const double scalar );
183 : : //- Compound Assignment: division: {this = this / scalar}
184 : :
185 : : VerdictVector operator-() const;
186 : : //- unary negation.
187 : :
188 : : friend VerdictVector operator~( const VerdictVector& vec );
189 : : //- normalize. Returns a new vector which is a copy of {vec},
190 : : //- scaled such that {|vec|=1}. Uses overloaded bitwise NOT operator.
191 : :
192 : : friend VerdictVector operator+( const VerdictVector& v1, const VerdictVector& v2 );
193 : : //- vector addition
194 : :
195 : : friend VerdictVector operator-( const VerdictVector& v1, const VerdictVector& v2 );
196 : : //- vector subtraction
197 : :
198 : : friend VerdictVector operator*( const VerdictVector& v1, const VerdictVector& v2 );
199 : : //- vector cross product, non-commutative
200 : :
201 : : friend VerdictVector operator*( const VerdictVector& v1, const double sclr );
202 : : //- vector * scalar
203 : :
204 : : friend VerdictVector operator*( const double sclr, const VerdictVector& v1 );
205 : : //- scalar * vector
206 : :
207 : : friend double operator%( const VerdictVector& v1, const VerdictVector& v2 );
208 : : //- dot product
209 : :
210 : : friend VerdictVector operator/( const VerdictVector& v1, const double sclr );
211 : : //- vector / scalar
212 : :
213 : : friend int operator==( const VerdictVector& v1, const VerdictVector& v2 );
214 : : //- Equality operator
215 : :
216 : : friend int operator!=( const VerdictVector& v1, const VerdictVector& v2 );
217 : : //- Inequality operator
218 : :
219 : : friend VerdictVector interpolate( const double param, const VerdictVector& v1, const VerdictVector& v2 );
220 : : //- Interpolate between two vectors. Returns (1-param)*v1 + param*v2
221 : :
222 : : VerdictVector& operator=( const VerdictVector& from );
223 : :
224 : : private:
225 : : double xVal; //- x component of vector.
226 : : double yVal; //- y component of vector.
227 : : double zVal; //- z component of vector.
228 : : };
229 : :
230 : : VerdictVector vectorRotate( const double angle, const VerdictVector& normalAxis, const VerdictVector& referenceAxis );
231 : : //- A new coordinate system is created with the xy plane corresponding
232 : : //- to the plane normal to {normalAxis}, and the x axis corresponding to
233 : : //- the projection of {referenceAxis} onto the normal plane. The normal
234 : : //- plane is the tangent plane at the root point. A unit vector is
235 : : //- constructed along the local x axis and then rotated by the given
236 : : //- ccw angle to form the new point. The new point, then is a unit
237 : : //- distance from the global origin in the tangent plane.
238 : : //- {angle} is in radians.
239 : :
240 : 49862 : inline double VerdictVector::x() const
241 : : {
242 : 49862 : return xVal;
243 : : }
244 : 49862 : inline double VerdictVector::y() const
245 : : {
246 : 49862 : return yVal;
247 : : }
248 : 49862 : inline double VerdictVector::z() const
249 : : {
250 : 49862 : return zVal;
251 : : }
252 : 132 : inline void VerdictVector::get_xyz( double xyz[3] )
253 : : {
254 : 132 : xyz[0] = xVal;
255 : 132 : xyz[1] = yVal;
256 : 132 : xyz[2] = zVal;
257 : 132 : }
258 : : inline void VerdictVector::get_xyz( double& xv, double& yv, double& zv )
259 : : {
260 : : xv = xVal;
261 : : yv = yVal;
262 : : zv = zVal;
263 : : }
264 : 0 : inline double& VerdictVector::r()
265 : : {
266 : 0 : return xVal;
267 : : }
268 : 0 : inline double& VerdictVector::theta()
269 : : {
270 : 0 : return yVal;
271 : : }
272 : 0 : inline void VerdictVector::x( const double xv )
273 : : {
274 : 0 : xVal = xv;
275 : 0 : }
276 : 0 : inline void VerdictVector::y( const double yv )
277 : : {
278 : 0 : yVal = yv;
279 : 0 : }
280 : 0 : inline void VerdictVector::z( const double zv )
281 : : {
282 : 0 : zVal = zv;
283 : 0 : }
284 : 0 : inline void VerdictVector::r( const double xv )
285 : : {
286 : 0 : xVal = xv;
287 : 0 : }
288 : 0 : inline void VerdictVector::theta( const double yv )
289 : : {
290 : 0 : yVal = yv;
291 : 0 : }
292 : 10936 : inline VerdictVector& VerdictVector::operator+=( const VerdictVector& vector )
293 : : {
294 : 10936 : xVal += vector.x();
295 : 10936 : yVal += vector.y();
296 : 10936 : zVal += vector.z();
297 : 10936 : return *this;
298 : : }
299 : :
300 : 1332 : inline VerdictVector& VerdictVector::operator-=( const VerdictVector& vector )
301 : : {
302 : 1332 : xVal -= vector.x();
303 : 1332 : yVal -= vector.y();
304 : 1332 : zVal -= vector.z();
305 : 1332 : return *this;
306 : : }
307 : :
308 : 7139 : inline VerdictVector& VerdictVector::operator*=( const VerdictVector& vector )
309 : : {
310 : : double xcross, ycross, zcross;
311 : 7139 : xcross = yVal * vector.z() - zVal * vector.y();
312 : 7139 : ycross = zVal * vector.x() - xVal * vector.z();
313 : 7139 : zcross = xVal * vector.y() - yVal * vector.x();
314 : 7139 : xVal = xcross;
315 : 7139 : yVal = ycross;
316 : 7139 : zVal = zcross;
317 : 7139 : return *this;
318 : : }
319 : :
320 : 35185 : inline VerdictVector::VerdictVector( const VerdictVector& copy_from )
321 : 35185 : : xVal( copy_from.xVal ), yVal( copy_from.yVal ), zVal( copy_from.zVal )
322 : : {
323 : 35185 : }
324 : :
325 : 8442 : inline VerdictVector::VerdictVector() : xVal( 0 ), yVal( 0 ), zVal( 0 ) {}
326 : :
327 : : inline VerdictVector::VerdictVector( const VerdictVector& tail, const VerdictVector& head )
328 : : : xVal( head.xVal - tail.xVal ), yVal( head.yVal - tail.yVal ), zVal( head.zVal - tail.zVal )
329 : : {
330 : : }
331 : :
332 : 4651 : inline VerdictVector::VerdictVector( const double xIn, const double yIn, const double zIn )
333 : 4651 : : xVal( xIn ), yVal( yIn ), zVal( zIn )
334 : : {
335 : 4651 : }
336 : :
337 : : // This sets the vector to be perpendicular to it's current direction.
338 : : // NOTE:
339 : : // This is a 2D function. It only works in the XY Plane.
340 : : inline void VerdictVector::perpendicular_z()
341 : : {
342 : : double temp = x();
343 : : x( y() );
344 : : y( -temp );
345 : : }
346 : :
347 : 11177 : inline void VerdictVector::set( const double xv, const double yv, const double zv )
348 : : {
349 : 11177 : xVal = xv;
350 : 11177 : yVal = yv;
351 : 11177 : zVal = zv;
352 : 11177 : }
353 : :
354 : 0 : inline void VerdictVector::set( const double xyz[3] )
355 : : {
356 : 0 : xVal = xyz[0];
357 : 0 : yVal = xyz[1];
358 : 0 : zVal = xyz[2];
359 : 0 : }
360 : :
361 : : inline void VerdictVector::set( const VerdictVector& tail, const VerdictVector& head )
362 : : {
363 : : xVal = head.xVal - tail.xVal;
364 : : yVal = head.yVal - tail.yVal;
365 : : zVal = head.zVal - tail.zVal;
366 : : }
367 : :
368 : 5191 : inline VerdictVector& VerdictVector::operator=( const VerdictVector& from )
369 : : {
370 : 5191 : xVal = from.xVal;
371 : 5191 : yVal = from.yVal;
372 : 5191 : zVal = from.zVal;
373 : 5191 : return *this;
374 : : }
375 : :
376 : : inline void VerdictVector::set( const VerdictVector& to_copy )
377 : : {
378 : : *this = to_copy;
379 : : }
380 : :
381 : : // Scale all values by scalar.
382 : 10534 : inline VerdictVector& VerdictVector::operator*=( const double scalar )
383 : : {
384 : 10534 : xVal *= scalar;
385 : 10534 : yVal *= scalar;
386 : 10534 : zVal *= scalar;
387 : 10534 : return *this;
388 : : }
389 : :
390 : : // Scales all values by 1/scalar
391 : 90 : inline VerdictVector& VerdictVector::operator/=( const double scalar )
392 : : {
393 [ - + ]: 90 : assert( scalar != 0 );
394 : 90 : xVal /= scalar;
395 : 90 : yVal /= scalar;
396 : 90 : zVal /= scalar;
397 : 90 : return *this;
398 : : }
399 : :
400 : : // Returns the normalized 'this'.
401 : : inline VerdictVector operator~( const VerdictVector& vec )
402 : : {
403 : : double mag = sqrt( vec.xVal * vec.xVal + vec.yVal * vec.yVal + vec.zVal * vec.zVal );
404 : :
405 : : VerdictVector temp = vec;
406 : : if( mag != 0.0 ) { temp /= mag; }
407 : : return temp;
408 : : }
409 : :
410 : : // Unary minus. Negates all values in vector.
411 : 445 : inline VerdictVector VerdictVector::operator-() const
412 : : {
413 : 445 : return VerdictVector( -xVal, -yVal, -zVal );
414 : : }
415 : :
416 : 330 : inline VerdictVector operator+( const VerdictVector& vector1, const VerdictVector& vector2 )
417 : : {
418 : 330 : double xv = vector1.x() + vector2.x();
419 : 330 : double yv = vector1.y() + vector2.y();
420 : 330 : double zv = vector1.z() + vector2.z();
421 : 330 : return VerdictVector( xv, yv, zv );
422 : : // return VerdictVector(vector1) += vector2;
423 : : }
424 : :
425 : 3346 : inline VerdictVector operator-( const VerdictVector& vector1, const VerdictVector& vector2 )
426 : : {
427 : 3346 : double xv = vector1.x() - vector2.x();
428 : 3346 : double yv = vector1.y() - vector2.y();
429 : 3346 : double zv = vector1.z() - vector2.z();
430 : 3346 : return VerdictVector( xv, yv, zv );
431 : : // return VerdictVector(vector1) -= vector2;
432 : : }
433 : :
434 : : // Cross products.
435 : : // vector1 cross vector2
436 : 7139 : inline VerdictVector operator*( const VerdictVector& vector1, const VerdictVector& vector2 )
437 : : {
438 [ + - ][ + - ]: 7139 : return VerdictVector( vector1 ) *= vector2;
439 : : }
440 : :
441 : : // Returns a scaled vector.
442 : : inline VerdictVector operator*( const VerdictVector& vector1, const double scalar )
443 : : {
444 : : return VerdictVector( vector1 ) *= scalar;
445 : : }
446 : :
447 : : // Returns a scaled vector
448 : 10274 : inline VerdictVector operator*( const double scalar, const VerdictVector& vector1 )
449 : : {
450 [ + - ][ + - ]: 10274 : return VerdictVector( vector1 ) *= scalar;
451 : : }
452 : :
453 : : // Returns a vector scaled by 1/scalar
454 : 90 : inline VerdictVector operator/( const VerdictVector& vector1, const double scalar )
455 : : {
456 [ + - ][ + - ]: 90 : return VerdictVector( vector1 ) /= scalar;
457 : : }
458 : :
459 : : inline int operator==( const VerdictVector& v1, const VerdictVector& v2 )
460 : : {
461 : : return ( v1.xVal == v2.xVal && v1.yVal == v2.yVal && v1.zVal == v2.zVal );
462 : : }
463 : :
464 : : inline int operator!=( const VerdictVector& v1, const VerdictVector& v2 )
465 : : {
466 : : return ( v1.xVal != v2.xVal || v1.yVal != v2.yVal || v1.zVal != v2.zVal );
467 : : }
468 : :
469 : 3479 : inline double VerdictVector::length_squared() const
470 : : {
471 : 3479 : return ( xVal * xVal + yVal * yVal + zVal * zVal );
472 : : }
473 : :
474 : 3181 : inline double VerdictVector::length() const
475 : : {
476 : 3181 : return ( sqrt( xVal * xVal + yVal * yVal + zVal * zVal ) );
477 : : }
478 : :
479 : 379 : inline double VerdictVector::normalize()
480 : : {
481 : 379 : double mag = length();
482 [ + - ]: 379 : if( mag != 0 )
483 : : {
484 : 379 : xVal = xVal / mag;
485 : 379 : yVal = yVal / mag;
486 : 379 : zVal = zVal / mag;
487 : : }
488 : 379 : return mag;
489 : : }
490 : :
491 : : // Dot Product.
492 : 7982 : inline double operator%( const VerdictVector& vector1, const VerdictVector& vector2 )
493 : : {
494 : 7982 : return ( vector1.x() * vector2.x() + vector1.y() * vector2.y() + vector1.z() * vector2.z() );
495 : : }
496 : :
497 : : #endif
|
