Line data Source code
1 : /*
2 : * Copyright (c) 2006-2015, Visualization and Multimedia Lab,
3 : * University of Zurich <http://vmml.ifi.uzh.ch>,
4 : * Eyescale Software GmbH,
5 : * Blue Brain Project, EPFL
6 : *
7 : * This file is part of VMMLib <https://github.com/VMML/vmmlib/>
8 : *
9 : * Redistribution and use in source and binary forms, with or without
10 : * modification, are permitted provided that the following conditions are met:
11 : *
12 : * Redistributions of source code must retain the above copyright notice, this
13 : * list of conditions and the following disclaimer. Redistributions in binary
14 : * form must reproduce the above copyright notice, this list of conditions and
15 : * the following disclaimer in the documentation and/or other materials provided
16 : * with the distribution. Neither the name of the Visualization and Multimedia
17 : * Lab, University of Zurich nor the names of its contributors may be used to
18 : * endorse or promote products derived from this software without specific prior
19 : * written permission.
20 : * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
21 : * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
22 : * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
23 : * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
24 : * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
25 : * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
26 : * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
27 : * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
28 : * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
29 : * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
30 : * POSSIBILITY OF SUCH DAMAGE.
31 : */
32 :
33 : #ifndef __VMML__VECTOR__HPP__
34 : #define __VMML__VECTOR__HPP__
35 :
36 : #include <vmmlib/vmmlib_config.hpp>
37 : #include <vmmlib/math.hpp>
38 : #include <vmmlib/enable_if.hpp>
39 : #include <vmmlib/exception.hpp>
40 :
41 : #include <iostream>
42 : #include <iomanip>
43 : #include <vector>
44 : #include <string>
45 : #include <cstring>
46 : #include <limits>
47 : #include <algorithm>
48 :
49 : namespace vmml
50 : {
51 :
52 : template< size_t M, typename T = float >
53 : class vector
54 : {
55 : public:
56 : typedef T value_type;
57 : typedef T* iterator;
58 : typedef const T* const_iterator;
59 : typedef std::reverse_iterator< iterator > reverse_iterator;
60 : typedef std::reverse_iterator< const_iterator > const_reverse_iterator;
61 :
62 : static const size_t DIMENSION = M;
63 :
64 : // constructors
65 79 : vector() : array() {} // http://stackoverflow.com/questions/5602030
66 : explicit vector( const T& a ); // sets all components to a;
67 : vector( const T& x, const T& y );
68 : vector( const T& x, const T& y, const T& z );
69 : vector( const T& x, const T& y, const T& z, const T& w );
70 :
71 : #ifndef SWIG
72 : // initializes the first M-1 values from vector_, the last from last_
73 : vector( const vector< M-1, T >& vector_, T last_ );
74 : #endif
75 :
76 : vector( const T* values );
77 :
78 : #ifdef __OSG_MATH
79 : template< typename OSGVEC3 >
80 : explicit vector( const OSGVEC3& from,
81 : typename enable_if< M == 3, OSGVEC3 >::type* = 0 );
82 : #endif
83 :
84 : // vec< M > with homogeneous coordinates <-> vec< M-1 > conversion ctor
85 : // to-homogenous-coordinates ctor
86 : template< size_t N >
87 : vector( const vector< N, T >& source_,
88 : typename enable_if< N == M - 1 >::type* = 0 );
89 :
90 : // from-homogenous-coordinates vector
91 : template< size_t N >
92 : vector( const vector< N, T >& source_,
93 : typename enable_if< N == M + 1 >::type* = 0 );
94 :
95 : template< typename U > vector( const vector< M, U >& source_ );
96 :
97 : // iterators
98 : inline iterator begin();
99 : inline iterator end();
100 : inline const_iterator begin() const;
101 : inline const_iterator end() const;
102 : inline reverse_iterator rbegin();
103 : inline reverse_iterator rend();
104 : inline const_reverse_iterator rbegin() const;
105 : inline const_reverse_iterator rend() const;
106 :
107 : # ifndef VMMLIB_NO_CONVERSION_OPERATORS
108 : // conversion operators
109 : inline operator T*();
110 : inline operator const T*() const;
111 : # else
112 : inline T& operator[]( size_t index );
113 : inline const T& operator[]( size_t index ) const;
114 : # endif
115 :
116 : // accessors
117 : inline T& operator()( size_t index );
118 : inline const T& operator()( size_t index ) const;
119 :
120 : inline T& at( size_t index );
121 : inline const T& at( size_t index ) const;
122 :
123 : // element accessors for M <= 4;
124 : inline T& x();
125 : inline T& y();
126 : inline T& z();
127 : inline T& w();
128 : inline const T& x() const;
129 : inline const T& y() const;
130 : inline const T& z() const;
131 : inline const T& w() const;
132 :
133 : // pixel color element accessors for M<= 4
134 : inline T& r();
135 : inline T& g();
136 : inline T& b();
137 : inline T& a();
138 : inline const T& r() const;
139 : inline const T& g() const;
140 : inline const T& b() const;
141 : inline const T& a() const;
142 :
143 : bool operator==( const vector& other ) const;
144 : bool operator!=( const vector& other ) const;
145 : bool equals( const vector& other,
146 : T tolerance = std::numeric_limits< T >::epsilon() ) const;
147 : bool operator<( const vector& other ) const;
148 :
149 : // remember kids: c_arrays are dangerous and evil!
150 : vector& operator=( const T* c_array );
151 : T operator=( T filler );
152 :
153 : vector& operator=( const vector& other );
154 :
155 : // returns void to avoid 'silent' loss of precision when chaining
156 : template< typename U > void operator=( const vector< M, U >& other );
157 :
158 : // to-homogenous-coordinates assignment operator
159 : // non-chainable because of sfinae
160 : template< size_t N >
161 : typename enable_if< N == M - 1 >::type*
162 : operator=( const vector< N, T >& source_ );
163 :
164 : // from-homogenous-coordinates assignment operator
165 : // non-chainable because of sfinae
166 : template< size_t N >
167 : typename enable_if< N == M + 1 >::type*
168 : operator=( const vector< N, T >& source_ );
169 :
170 : vector operator*( const vector& other ) const;
171 : vector operator/( const vector& other ) const;
172 : vector operator+( const vector& other ) const;
173 : vector operator-( const vector& other ) const;
174 :
175 : void operator*=( const vector& other );
176 : void operator/=( const vector& other );
177 : void operator+=( const vector& other );
178 : void operator-=( const vector& other );
179 :
180 : vector operator*( const T other ) const;
181 : vector operator/( const T other ) const;
182 : vector operator+( const T other ) const;
183 : vector operator-( const T other ) const;
184 :
185 : void operator*=( const T other );
186 : void operator/=( const T other );
187 : void operator+=( const T other );
188 : void operator-=( const T other );
189 :
190 : vector operator-() const;
191 :
192 : const vector& negate();
193 :
194 : void set( T a ); // sets all components to a;
195 : #ifndef SWIG
196 : void set( const vector< M-1, T >& v, T a );
197 : #endif
198 : template< size_t N >
199 : void set( const vector< N, T >& v );
200 :
201 : // sets the first few components to a certain value
202 : void set( T x, T y );
203 : void set( T x, T y, T z );
204 : void set( T x, T y, T z, T w );
205 :
206 : template< typename input_iterator_t >
207 : void iter_set( input_iterator_t begin_, input_iterator_t end_ );
208 :
209 : // compute the cross product of two vectors
210 : // note: there's also a free function:
211 : // vector<> cross( const vector<>, const vector<> )
212 :
213 : // result = vec1.cross( vec2 ) => retval result = vec1 x vec2
214 : template< typename TT >
215 : vector cross( const vector< M, TT >& rhs,
216 : typename enable_if< M == 3, TT >::type* = 0 ) const;
217 :
218 : // result.cross( vec1, vec2 ) => (this) = vec1 x vec2
219 : template< typename TT >
220 : void cross( const vector< M, TT >& a, const vector< M, TT >& b,
221 : typename enable_if< M == 3, TT >::type* = 0 );
222 :
223 :
224 : // compute the dot product of two vectors
225 : // note: there's also a free function:
226 : // T dot( const vector<>, const vector<> );
227 : inline T dot( const vector& other ) const;
228 :
229 :
230 : // normalize the vector
231 : // note: there's also a free function:
232 : // vector<> normalize( const vector<> );
233 : inline T normalize();
234 :
235 : //sets all vector components to random values
236 : //remember to set srand( seed );
237 : //if seed is set to -1, srand( seed ) was set outside set_random
238 : //otherwise srand( seed ) will be called with the given seed
239 : void set_random( int seed = -1 );
240 :
241 : inline T length() const;
242 : inline T squared_length() const;
243 :
244 : inline T distance( const vector& other ) const;
245 : inline T squared_distance( const vector& other ) const;
246 :
247 : /** @return the product of all elements of this vector */
248 : T product() const;
249 :
250 : template< typename TT >
251 : vector< 3, T > rotate( const T theta, vector< M, TT > axis,
252 : typename enable_if< M == 3, TT >::type* = 0 ) const;
253 :
254 : // right hand system, CCW triangle
255 : // (*this) = normal of v0,v1,v2
256 : void compute_normal( const vector& v0, const vector& v1, const vector& v2 );
257 : // retval = normal of (this), v1, v2
258 : vector compute_normal( const vector& v1, const vector& v2 ) const;
259 :
260 : /** @return the sub vector at the given position and length. */
261 : template< size_t N >
262 : vector< N, T >& get_sub_vector( size_t offset = 0,
263 : typename enable_if< M >= N >::type* = 0 );
264 :
265 : /** @return the sub vector at the given position and length. */
266 : template< size_t N >
267 : const vector< N, T >& get_sub_vector( size_t offset = 0,
268 : typename enable_if< M >= N >::type* = 0 ) const;
269 :
270 : // sphere functions - sphere layout: center xyz, radius w
271 : template< typename TT >
272 : inline vector< 3, T > project_point_onto_sphere(
273 : const vector< 3, TT >& point,
274 : typename enable_if< M == 4, TT >::type* = 0 ) const;
275 :
276 : // returns a negative distance if the point lies in the sphere
277 : template< typename TT >
278 : inline T distance_to_sphere( const vector< 3, TT >& point,
279 : typename enable_if< M == 4, TT >::type* = 0 ) const;
280 :
281 : // plane functions - plane layout; normal xyz, distance w
282 : template< typename TT >
283 : inline T distance_to_plane( const vector< 3, TT >& point,
284 : typename enable_if< M == 4, TT >::type* = 0 ) const;
285 :
286 : template< typename TT >
287 : inline vector< 3, T > project_point_onto_plane(
288 : const vector< 3, TT >& point,
289 : typename enable_if< M == 4, TT >::type* = 0 ) const;
290 :
291 : // returns the index of the minimal resp. maximal value in the vector
292 : size_t find_min_index() const;
293 : size_t find_max_index() const;
294 :
295 : // returns the index of the minimal resp. maximal value in the vector
296 : size_t find_abs_min_index() const;
297 : size_t find_abs_max_index() const;
298 :
299 : // returns minimal resp. maximal value in the vector
300 : T& find_min();
301 : T& find_max();
302 : const T& find_min() const;
303 : const T& find_max() const;
304 :
305 : void clamp( const T& min = 0.0, const T& max = 1.0 );
306 :
307 : template< typename TT >
308 : void scale_to( vector< M, TT >& scaled_vector,
309 : T min_value = -1.0, T max_value = 1.0 ) const;
310 :
311 : inline static size_t size(); // returns M
312 :
313 : bool is_unit_vector() const;
314 :
315 : // perturbs each component by randomly + or - the perturbation parameter
316 : void perturb( T perturbation = 0.0001 );
317 :
318 : void sqrt_elementwise();
319 : double norm() const; //l2 norm
320 :
321 : // computes the reciprocal value for each component, x = 1/x;
322 : // WARNING: might result in nans if division by 0!
323 : void reciprocal();
324 : // computes the reciprocal value for each component, x = 1/x;
325 : // checks every component for 0, sets to max value if zero.
326 : void reciprocal_safe();
327 :
328 : template< typename TT >
329 : void cast_from( const vector< M, TT >& other );
330 :
331 : size_t nnz() const;
332 :
333 : // test each component of the vector for isnan and isinf
334 : // inline bool is_valid() const; -> moved to class validator
335 :
336 0 : friend std::ostream& operator<< ( std::ostream& os, const vector& vector_ )
337 : {
338 0 : const std::ios::fmtflags flags = os.flags();
339 0 : const int prec = os.precision();
340 :
341 0 : os.setf( std::ios::right, std::ios::adjustfield );
342 0 : os.precision( 5 );
343 0 : os << "[ ";
344 0 : for( size_t index = 0; index < M; ++index )
345 0 : os << std::setw(10) << vector_.at( index ) << " ";
346 0 : os << "]";
347 0 : os.precision( prec );
348 0 : os.setf( flags );
349 0 : return os;
350 : }
351 :
352 : T array[ M ]; //!< storage
353 :
354 : #ifndef SWIG
355 : // Vector3 defaults
356 : static const vector FORWARD;
357 : static const vector BACKWARD;
358 : static const vector UP;
359 : static const vector DOWN;
360 : static const vector LEFT;
361 : static const vector RIGHT;
362 :
363 : static const vector ONE;
364 : static const vector ZERO;
365 :
366 : // Unit vectors
367 : static const vector UNIT_X;
368 : static const vector UNIT_Y;
369 : static const vector UNIT_Z;
370 : #endif
371 :
372 : }; // class vector
373 :
374 : //
375 : // typedefs and statics
376 : //
377 : #ifndef SWIG
378 : template< size_t M, typename T >
379 : const vector< M, T > vector< M, T >::FORWARD( 0, 0, -1 );
380 : template< size_t M, typename T >
381 : const vector< M, T > vector< M, T >::BACKWARD( 0, 0, 1 );
382 : template< size_t M, typename T >
383 : const vector< M, T > vector< M, T >::UP( 0, 1, 0 );
384 : template< size_t M, typename T >
385 : const vector< M, T > vector< M, T >::DOWN( 0, -1, 0 );
386 : template< size_t M, typename T >
387 : const vector< M, T > vector< M, T >::LEFT( -1, 0, 0 );
388 : template< size_t M, typename T >
389 : const vector< M, T > vector< M, T >::RIGHT( 1, 0, 0 );
390 : template< size_t M, typename T >
391 :
392 : const vector< M, T > vector< M, T >::ONE( static_cast< T >( 1 ));
393 : template< size_t M, typename T >
394 : const vector< M, T > vector< M, T >::ZERO( static_cast< T >( 0 ));
395 : template< size_t M, typename T >
396 :
397 : const vector< M, T > vector< M, T >::UNIT_X( 1, 0, 0 );
398 : template< size_t M, typename T >
399 : const vector< M, T > vector< M, T >::UNIT_Y( 0, 1, 0 );
400 : template< size_t M, typename T >
401 : const vector< M, T > vector< M, T >::UNIT_Z( 0, 0, 1 );
402 : #endif
403 :
404 : #ifndef VMMLIB_NO_TYPEDEFS
405 : # ifdef _MSC_VER
406 : typedef UINT8 uint8_t;
407 : # endif
408 : typedef vmml::vector< 2, int > Vector2i;
409 : typedef vmml::vector< 3, int > Vector3i;
410 : typedef vmml::vector< 4, int > Vector4i;
411 : typedef vmml::vector< 2, unsigned > Vector2ui;
412 : typedef vmml::vector< 3, unsigned > Vector3ui;
413 : typedef vmml::vector< 4, unsigned > Vector4ui;
414 : typedef vmml::vector< 3, double > Vector3d;
415 : typedef vmml::vector< 4, double > Vector4d;
416 : typedef vmml::vector< 2, float > Vector2f;
417 : typedef vmml::vector< 3, float > Vector3f;
418 : typedef vmml::vector< 4, float > Vector4f;
419 : typedef vmml::vector< 3, uint8_t > Vector3ub;
420 : typedef vmml::vector< 4, uint8_t > Vector4ub;
421 : #endif
422 :
423 : //
424 : // some free functions for convenience
425 : //
426 :
427 : template< size_t M, typename T >
428 : bool equals( const vector< M, T >& a, const vector< M, T >& b )
429 : {
430 : return a.equals( b );
431 : }
432 :
433 :
434 : // allows float * vector, not only vector * float
435 : template< size_t M, typename T >
436 1 : static vector< M, T > operator* ( T factor, const vector< M, T >& vector_ )
437 : {
438 1 : return vector_ * factor;
439 : }
440 :
441 :
442 : template< size_t M, typename T >
443 : inline T dot( const vector< M, T >& first, const vector< M, T >& second )
444 : {
445 : return first.dot( second );
446 : }
447 :
448 :
449 : template< size_t M, typename T >
450 : inline vector< M, T > cross( const vector< 3, T >& a, const vector< 3, T >& b )
451 : {
452 : return a.cross( b );
453 : }
454 :
455 :
456 : template< size_t M, typename T >
457 : inline vector< M, T > normalize( const vector< M, T >& vector_ )
458 : {
459 : vector< M, T > v( vector_ );
460 : v.normalize();
461 : return v;
462 : }
463 :
464 : template< typename T >
465 6 : inline vector< 4, T > compute_plane( const vector< 3, T >& a,
466 : const vector< 3, T >& b,
467 : const vector< 3, T >& c )
468 : {
469 6 : const vector< 3, T > cb = b - c;
470 6 : const vector< 3, T > ba = a - b;
471 :
472 6 : vector< 4, T > plane = cb.cross( ba );
473 6 : plane.normalize();
474 6 : plane.w() = -plane.x() * a.x() - plane.y() * a.y() - plane.z() * a.z();
475 6 : return plane;
476 : }
477 :
478 : template< size_t M, typename T >
479 3 : vector< M, T >::vector( const T& _a )
480 : {
481 13 : for( iterator it = begin(), it_end = end(); it != it_end; ++it )
482 : {
483 10 : *it = _a;
484 : }
485 3 : }
486 :
487 : template< size_t M, typename T >
488 9 : vector< M, T >::vector( const T& _x, const T& _y )
489 : {
490 9 : array[ 0 ] = _x;
491 9 : array[ 1 ] = _y;
492 9 : }
493 :
494 :
495 : template< size_t M, typename T >
496 33 : vector< M, T >::vector( const T& _x, const T& _y, const T& _z )
497 : {
498 33 : array[ 0 ] = _x;
499 33 : array[ 1 ] = _y;
500 33 : array[ 2 ] = _z;
501 33 : }
502 :
503 :
504 :
505 : template< size_t M, typename T >
506 19 : vector< M, T >::vector( const T& _x, const T& _y, const T& _z, const T& _w )
507 : {
508 19 : array[ 0 ] = _x;
509 19 : array[ 1 ] = _y;
510 19 : array[ 2 ] = _z;
511 19 : array[ 3 ] = _w;
512 19 : }
513 :
514 :
515 : template< size_t M, typename T >
516 : vector< M, T >::vector( const T* values )
517 : {
518 : memcpy( array, values, M * sizeof( T ));
519 : }
520 :
521 : #ifdef __OSG_MATH
522 : template< size_t M, typename T >
523 : template< typename OSGVEC3 >
524 : vector< M, T >::vector( const OSGVEC3& from,
525 : typename enable_if< M == 3, OSGVEC3 >::type* )
526 : {
527 : array[ 0 ] = from.x();
528 : array[ 1 ] = from.y();
529 : array[ 2 ] = from.z();
530 : }
531 : #endif
532 :
533 : #ifndef SWIG
534 : template< size_t M, typename T >
535 : // initializes the first M-1 values from vector_, the last from last_
536 1 : vector< M, T >::vector( const vector< M-1, T >& vector_, T last_ )
537 : {
538 : typename vector< M-1, T >::const_iterator
539 1 : it = vector_.begin(), it_end = vector_.end();
540 :
541 1 : iterator my_it = begin();
542 :
543 4 : for( ; it != it_end; ++it, ++my_it )
544 : {
545 3 : (*my_it) = *it;
546 : }
547 1 : (*my_it) = last_;
548 1 : }
549 : #endif
550 :
551 :
552 :
553 : // to-homogenous-coordinates ctor
554 : template< size_t M, typename T >
555 : template< size_t N >
556 41 : vector< M, T >::vector( const vector< N, T >& source_,
557 : typename enable_if< N == M - 1 >::type* )
558 : {
559 41 : (*this) = source_;
560 41 : }
561 :
562 :
563 :
564 :
565 : // from-homogenous-coordinates ctor
566 : template< size_t M, typename T >
567 : template< size_t N >
568 1 : vector< M, T >::vector( const vector< N, T >& source_,
569 : typename enable_if< N == M + 1 >::type* )
570 : {
571 1 : (*this) = source_;
572 1 : }
573 :
574 :
575 : template< size_t M, typename T >
576 : template< typename U >
577 1 : vector< M, T >::vector( const vector< M, U >& source_ )
578 : {
579 1 : (*this) = source_;
580 1 : }
581 :
582 :
583 :
584 1 : template< size_t M, typename T > void vector< M, T >::set( T _a )
585 : {
586 5 : for( iterator it = begin(), it_end = end(); it != it_end; ++it )
587 4 : *it = _a;
588 1 : }
589 :
590 :
591 : #ifndef SWIG
592 : template< size_t M, typename T >
593 : void vector< M, T >::set( const vector< M-1, T >& v, T _a )
594 : {
595 : memcpy( array, v.array, sizeof( T ) * (M-1) );
596 : at( M-1 ) = _a;
597 : }
598 : #endif
599 :
600 : template< size_t M, typename T > template< size_t N >
601 : void vector< M, T >::set( const vector< N, T >& v )
602 : {
603 : size_t minimum = M;
604 : if (N < M) minimum = N;
605 : memcpy( array, v.array, sizeof( T ) * minimum );
606 : }
607 :
608 : template< size_t M, typename T >
609 : void vector< M, T >::set( T _x, T _y )
610 : {
611 : array[ 0 ] = _x;
612 : array[ 1 ] = _y;
613 : }
614 :
615 :
616 : template< size_t M, typename T >
617 : void vector< M, T >::set( T _x, T _y, T _z )
618 : {
619 : array[ 0 ] = _x;
620 : array[ 1 ] = _y;
621 : array[ 2 ] = _z;
622 : }
623 :
624 :
625 :
626 : template< size_t M, typename T >
627 1 : void vector< M, T >::set( T _x, T _y, T _z, T _w )
628 : {
629 1 : array[ 0 ] = _x;
630 1 : array[ 1 ] = _y;
631 1 : array[ 2 ] = _z;
632 1 : array[ 3 ] = _w;
633 1 : }
634 :
635 :
636 : template< size_t M, typename T >
637 : inline T&
638 : vector< M, T >::operator()( size_t index )
639 : {
640 : return at( index );
641 : }
642 :
643 :
644 :
645 : template< size_t M, typename T >
646 : inline const T&
647 4 : vector< M, T >::operator()( size_t index ) const
648 : {
649 4 : return at( index );
650 : }
651 :
652 :
653 :
654 : template< size_t M, typename T >
655 : inline T&
656 346 : vector< M, T >::at( size_t index )
657 : {
658 : #ifdef VMMLIB_SAFE_ACCESSORS
659 346 : if ( index >= M )
660 : {
661 0 : VMMLIB_ERROR( "at() - index out of bounds", VMMLIB_HERE );
662 : }
663 : #endif
664 346 : return array[ index ];
665 : }
666 :
667 :
668 :
669 : template< size_t M, typename T >
670 : inline const T&
671 524 : vector< M, T >::at( size_t index ) const
672 : {
673 : #ifdef VMMLIB_SAFE_ACCESSORS
674 524 : if ( index >= M )
675 : {
676 0 : VMMLIB_ERROR( "at() - index out of bounds", VMMLIB_HERE );
677 : }
678 : #endif
679 524 : return array[ index ];
680 : }
681 :
682 :
683 : #ifndef VMMLIB_NO_CONVERSION_OPERATORS
684 :
685 : template< size_t M, typename T >
686 : vector< M, T >::operator T*()
687 : {
688 : return array;
689 : }
690 :
691 :
692 :
693 : template< size_t M, typename T >
694 72 : vector< M, T >::operator const T*() const
695 : {
696 72 : return array;
697 : }
698 : #else
699 :
700 : template< size_t M, typename T >
701 : T&
702 : vector< M, T >::operator[]( size_t index )
703 : {
704 : return at( index );
705 : }
706 :
707 : template< size_t M, typename T >
708 : const T&
709 : vector< M, T >::operator[]( size_t index ) const
710 : {
711 : return at( index );
712 : }
713 :
714 :
715 : #endif
716 :
717 :
718 : #if 0
719 : template< size_t M, typename T >
720 : inline T&
721 : vector< M, T >::operator[]( size_t index )
722 : {
723 : return at( index );
724 : }
725 :
726 :
727 :
728 : template< size_t M, typename T >
729 : inline const T&
730 : vector< M, T >::operator[]( size_t index ) const
731 : {
732 : return at( index );
733 : }
734 : #endif
735 :
736 :
737 : template< size_t M, typename T >
738 : vector< M, T >
739 1 : vector< M, T >::operator*( const vector< M, T >& other ) const
740 : {
741 1 : vector< M, T > result;
742 5 : for( size_t index = 0; index < M; ++index )
743 4 : result.at( index ) = at( index ) * other.at( index );
744 1 : return result;
745 : }
746 :
747 :
748 :
749 : template< size_t M, typename T >
750 : vector< M, T >
751 1 : vector< M, T >::operator/( const vector< M, T >& other ) const
752 : {
753 1 : vector< M, T > result;
754 5 : for( size_t index = 0; index < M; ++index )
755 4 : result.at( index ) = at( index ) / other.at( index );
756 1 : return result;
757 : }
758 :
759 :
760 :
761 : template< size_t M, typename T >
762 : vector< M, T >
763 4 : vector< M, T >::operator+( const vector< M, T >& other ) const
764 : {
765 4 : vector< M, T > result;
766 20 : for( size_t index = 0; index < M; ++index )
767 16 : result.at( index ) = at( index ) + other.at( index );
768 4 : return result;
769 : }
770 :
771 :
772 :
773 : template< size_t M, typename T >
774 : vector< M, T >
775 18 : vector< M, T >::operator-( const vector< M, T >& other ) const
776 : {
777 18 : vector< M, T > result;
778 76 : for( size_t index = 0; index < M; ++index )
779 58 : result.at( index ) = at( index ) - other.at( index );
780 18 : return result;
781 : }
782 :
783 :
784 :
785 :
786 : template< size_t M, typename T >
787 : void
788 1 : vector< M, T >::operator*=( const vector< M, T >& other )
789 : {
790 5 : for( size_t index = 0; index < M; ++index )
791 4 : at( index ) *= other.at( index );
792 1 : }
793 :
794 :
795 :
796 : template< size_t M, typename T >
797 : void
798 1 : vector< M, T >::operator/=( const vector< M, T >& other )
799 : {
800 5 : for( size_t index = 0; index < M; ++index )
801 4 : at( index ) /= other.at( index );
802 1 : }
803 :
804 :
805 :
806 : template< size_t M, typename T >
807 : void
808 1 : vector< M, T >::operator+=( const vector< M, T >& other )
809 : {
810 5 : for( size_t index = 0; index < M; ++index )
811 4 : at( index ) += other.at( index );
812 1 : }
813 :
814 :
815 :
816 : template< size_t M, typename T >
817 : void
818 1 : vector< M, T >::operator-=( const vector< M, T >& other )
819 : {
820 5 : for( size_t index = 0; index < M; ++index )
821 4 : at( index ) -= other.at( index );
822 1 : }
823 :
824 :
825 :
826 : template< size_t M, typename T >
827 : vector< M, T >
828 14 : vector< M, T >::operator*( const T other ) const
829 : {
830 14 : vector< M, T > result;
831 58 : for( size_t index = 0; index < M; ++index )
832 44 : result.at( index ) = at( index ) * other;
833 14 : return result;
834 : }
835 :
836 :
837 :
838 : template< size_t M, typename T >
839 : vector< M, T >
840 1 : vector< M, T >::operator/( const T other ) const
841 : {
842 1 : vector< M, T > result;
843 5 : for( size_t index = 0; index < M; ++index )
844 4 : result.at( index ) = at( index ) / other;
845 1 : return result;
846 : }
847 :
848 :
849 :
850 : template< size_t M, typename T >
851 : vector< M, T >
852 1 : vector< M, T >::operator+( const T other ) const
853 : {
854 1 : vector< M, T > result;
855 5 : for( size_t index = 0; index < M; ++index )
856 4 : result.at( index ) = at( index ) + other;
857 1 : return result;
858 : }
859 :
860 :
861 :
862 : template< size_t M, typename T >
863 : vector< M, T >
864 1 : vector< M, T >::operator-( const T other ) const
865 : {
866 1 : vector< M, T > result;
867 5 : for( size_t index = 0; index < M; ++index )
868 4 : result.at( index ) = at( index ) - other;
869 1 : return result;
870 : }
871 :
872 :
873 :
874 :
875 : template< size_t M, typename T >
876 : void
877 9 : vector< M, T >::operator*=( const T other )
878 : {
879 44 : for( size_t index = 0; index < M; ++index )
880 35 : at( index ) *= other;
881 9 : }
882 :
883 :
884 :
885 : template< size_t M, typename T >
886 : void
887 1 : vector< M, T >::operator/=( const T other )
888 : {
889 5 : for( size_t index = 0; index < M; ++index )
890 4 : at( index ) /= other;
891 1 : }
892 :
893 :
894 :
895 : template< size_t M, typename T >
896 : void
897 1 : vector< M, T >::operator+=( const T other )
898 : {
899 5 : for( size_t index = 0; index < M; ++index )
900 4 : at( index ) += other;
901 1 : }
902 :
903 :
904 :
905 : template< size_t M, typename T >
906 : void
907 1 : vector< M, T >::operator-=( const T other )
908 : {
909 5 : for( size_t index = 0; index < M; ++index )
910 4 : at( index ) -= other;
911 1 : }
912 :
913 :
914 :
915 : template< size_t M, typename T >
916 : vector< M, T >
917 2 : vector< M, T >::operator-() const
918 : {
919 2 : vector< M, T > v( *this );
920 2 : return v.negate();
921 : }
922 :
923 :
924 :
925 : template< size_t M, typename T >
926 : const vector< M, T >&
927 2 : vector< M, T >::negate()
928 : {
929 8 : for( size_t index = 0; index < M; ++index )
930 6 : array[ index ] = -array[ index ];
931 2 : return *this;
932 : }
933 :
934 :
935 :
936 : template< size_t M, typename T >
937 : inline T&
938 27 : vector< M, T >::x()
939 : {
940 27 : return array[ 0 ];
941 : }
942 :
943 :
944 :
945 : template< size_t M, typename T >
946 : inline T&
947 27 : vector< M, T >::y()
948 : {
949 27 : return array[ 1 ];
950 : }
951 :
952 :
953 :
954 : template< size_t M, typename T >
955 : inline T&
956 23 : vector< M, T >::z()
957 : {
958 23 : return array[ 2 ];
959 : }
960 :
961 :
962 :
963 : template< size_t M, typename T >
964 : inline T&
965 13 : vector< M, T >::w()
966 : {
967 13 : return array[ 3 ];
968 : }
969 :
970 :
971 :
972 : template< size_t M, typename T >
973 : inline const T&
974 190 : vector< M, T >::x() const
975 : {
976 190 : return array[ 0 ];
977 : }
978 :
979 :
980 :
981 : template< size_t M, typename T >
982 : inline const T&
983 184 : vector< M, T >::y() const
984 : {
985 184 : return array[ 1 ];
986 : }
987 :
988 :
989 :
990 : template< size_t M, typename T >
991 : inline const T&
992 182 : vector< M, T >::z() const
993 : {
994 182 : return array[ 2 ];
995 : }
996 :
997 :
998 :
999 : template< size_t M, typename T >
1000 : inline const T&
1001 100 : vector< M, T >::w() const
1002 : {
1003 100 : return array[ 3 ];
1004 : }
1005 :
1006 :
1007 : template< size_t M, typename T >
1008 : inline T&
1009 : vector< M, T >::r()
1010 : {
1011 : return array[ 0 ];
1012 : }
1013 :
1014 :
1015 :
1016 : template< size_t M, typename T >
1017 : inline T&
1018 : vector< M, T >::g()
1019 : {
1020 : return array[ 1 ];
1021 : }
1022 :
1023 :
1024 :
1025 : template< size_t M, typename T >
1026 : inline T&
1027 : vector< M, T >::b()
1028 : {
1029 : return array[ 2 ];
1030 : }
1031 :
1032 :
1033 :
1034 : template< size_t M, typename T >
1035 : inline T&
1036 : vector< M, T >::a()
1037 : {
1038 : return array[ 3 ];
1039 : }
1040 :
1041 :
1042 :
1043 : template< size_t M, typename T >
1044 : inline const T&
1045 : vector< M, T >::r() const
1046 : {
1047 : return array[ 0 ];
1048 : }
1049 :
1050 :
1051 :
1052 : template< size_t M, typename T >
1053 : inline const T&
1054 : vector< M, T >::g() const
1055 : {
1056 : return array[ 1 ];
1057 : }
1058 :
1059 :
1060 :
1061 : template< size_t M, typename T >
1062 : inline const T&
1063 : vector< M, T >::b() const
1064 : {
1065 : return array[ 2 ];
1066 : }
1067 :
1068 :
1069 :
1070 : template< size_t M, typename T >
1071 : inline const T&
1072 : vector< M, T >::a() const
1073 : {
1074 : return array[ 3 ];
1075 : }
1076 :
1077 : // result = vec1.cross( vec2 ) => result = vec1 x vec2
1078 : template< size_t M, typename T >
1079 : template< typename TT >
1080 7 : inline vector< M, T > vector< M, T >::cross( const vector< M, TT >& rhs,
1081 : typename enable_if< M == 3, TT >::type* ) const
1082 : {
1083 7 : vector< M, T > result;
1084 7 : result.cross( *this, rhs );
1085 7 : return result;
1086 : }
1087 :
1088 :
1089 :
1090 : // result.cross( vec1, vec2 ) => (this) = vec1 x vec2
1091 : template< size_t M, typename T >
1092 : template< typename TT >
1093 7 : void vector< M, T >::cross( const vector< M, TT >& aa,
1094 : const vector< M, TT >& bb,
1095 : typename enable_if< M == 3, TT >::type* )
1096 : {
1097 7 : array[ 0 ] = aa.y() * bb.z() - aa.z() * bb.y();
1098 7 : array[ 1 ] = aa.z() * bb.x() - aa.x() * bb.z();
1099 7 : array[ 2 ] = aa.x() * bb.y() - aa.y() * bb.x();
1100 7 : }
1101 :
1102 :
1103 :
1104 : template< size_t M, typename T >
1105 35 : inline T vector< M, T >::dot( const vector< M, T >& other ) const
1106 : {
1107 35 : T tmp = 0.0;
1108 174 : for( size_t index = 0; index < M; ++index )
1109 139 : tmp += at( index ) * other.at( index );
1110 :
1111 35 : return tmp;
1112 : }
1113 :
1114 :
1115 : template< size_t M, typename T >
1116 8 : inline T vector< M, T >::normalize()
1117 : {
1118 8 : const T len = length();
1119 :
1120 8 : if ( len == 0 )
1121 0 : return 0;
1122 :
1123 8 : const T tmp = 1.0 / len;
1124 8 : (*this) *= tmp;
1125 8 : return len;
1126 : }
1127 :
1128 : template< size_t M, typename T >
1129 16 : inline T vector< M, T >::length() const
1130 : {
1131 16 : return std::sqrt( squared_length() );
1132 : }
1133 :
1134 : template< size_t M, typename T >
1135 17 : inline T vector< M, T >::squared_length() const
1136 : {
1137 17 : T _squared_length = 0.0;
1138 78 : for( const_iterator it = begin(), it_end = end(); it != it_end; ++it )
1139 61 : _squared_length += (*it) * (*it);
1140 :
1141 17 : return _squared_length;
1142 : }
1143 :
1144 :
1145 :
1146 : template< size_t M, typename T >
1147 : inline T
1148 : vector< M, T >::distance( const vector< M, T >& other ) const
1149 : {
1150 : return std::sqrt( squared_distance( other ) );
1151 : }
1152 :
1153 :
1154 :
1155 : template< size_t M, typename T >
1156 : inline T vector< M, T >::squared_distance( const vector< M, T >& other ) const
1157 : {
1158 : vector< M, T > tmp( *this );
1159 : tmp -= other;
1160 : return tmp.squared_length();
1161 : }
1162 :
1163 2 : template< size_t M, typename T > inline T vector< M, T >::product() const
1164 : {
1165 2 : T result = at( 0 );
1166 6 : for( size_t i = 1; i < M; ++i )
1167 4 : result *= at( i );
1168 2 : return result;
1169 : }
1170 :
1171 : template< size_t M, typename T >
1172 : void vector< M, T >::compute_normal( const vector< M, T >& aa,
1173 : const vector< M, T >& bb,
1174 : const vector< M, T >& cc )
1175 : {
1176 : vector< M, T > u,v;
1177 : // right hand system, CCW triangle
1178 : u = bb - aa;
1179 : v = cc - aa;
1180 : cross( u, v );
1181 : normalize();
1182 : }
1183 :
1184 :
1185 :
1186 : template< size_t M, typename T >
1187 : vector< M, T >
1188 : vector< M, T >::compute_normal( const vector< M, T >& bb,
1189 : const vector< M, T >& cc ) const
1190 : {
1191 : vector< M, T > tmp;
1192 : tmp.compute_normal( *this, bb, cc);
1193 : return tmp;
1194 : }
1195 :
1196 : template< size_t M, typename T >
1197 : template< typename TT >
1198 : vector< 3, T > vector< M, T >::rotate( const T theta, vector< M, TT > axis,
1199 : typename enable_if< M == 3, TT >::type* ) const
1200 : {
1201 : axis.normalize();
1202 :
1203 : const T costheta = std::cos( theta );
1204 : const T sintheta = std::sin( theta );
1205 :
1206 : return vector< 3, T >(
1207 : (costheta + ( 1.0f - costheta ) * axis.x() * axis.x() ) * x() +
1208 : (( 1 - costheta ) * axis.x() * axis.y() - axis.z() * sintheta ) * y() +
1209 : (( 1 - costheta ) * axis.x() * axis.z() + axis.y() * sintheta ) * z(),
1210 :
1211 : (( 1 - costheta ) * axis.x() * axis.y() + axis.z() * sintheta ) * x() +
1212 : ( costheta + ( 1 - costheta ) * axis.y() * axis.y() ) * y() +
1213 : (( 1 - costheta ) * axis.y() * axis.z() - axis.x() * sintheta ) * z(),
1214 :
1215 : (( 1 - costheta ) * axis.x() * axis.z() - axis.y() * sintheta ) * x() +
1216 : (( 1 - costheta ) * axis.y() * axis.z() + axis.x() * sintheta ) * y() +
1217 : ( costheta + ( 1 - costheta ) * axis.z() * axis.z() ) * z() );
1218 : }
1219 :
1220 :
1221 : // sphere layout: center xyz, radius w
1222 : template< size_t M, typename T >
1223 : template< typename TT >
1224 : inline vector< 3, T >
1225 : vector< M, T >::
1226 : project_point_onto_sphere( const vector< 3, TT >& point,
1227 : typename enable_if< M == 4, TT >::type* ) const
1228 : {
1229 : const vector< 3, T >& _center = get_sub_vector< 3 >( 0 );
1230 :
1231 : vector< 3, T > projected_point( point );
1232 : projected_point -= _center;
1233 : projected_point.normalize();
1234 : projected_point *= w();
1235 : return _center + projected_point;
1236 : }
1237 :
1238 :
1239 :
1240 : // sphere layout: center xyz, radius w
1241 : template< size_t M, typename T >
1242 : template< typename TT >
1243 : inline T
1244 : vector< M, T >::
1245 : distance_to_sphere( const vector< 3, TT >& point,
1246 : typename enable_if< M == 4, TT >::type* ) const
1247 : {
1248 : const vector< 3, T >& center_ = get_sub_vector< 3 >( 0 );
1249 : return ( point - center_ ).length() - w();
1250 : }
1251 :
1252 : template< size_t M, typename T > template< size_t N > inline
1253 7 : vector< N, T >& vector< M, T >::get_sub_vector( size_t offset,
1254 : typename enable_if< M >= N >::type* )
1255 : {
1256 7 : assert( offset <= M - N );
1257 7 : return reinterpret_cast< vector< N, T >& >( *( begin() + offset ) );
1258 : }
1259 :
1260 : template< size_t M, typename T > template< size_t N > inline
1261 : const vector< N, T >& vector< M, T >::get_sub_vector( size_t offset,
1262 : typename enable_if< M >= N >::type* ) const
1263 : {
1264 : assert( offset <= M - N );
1265 : return reinterpret_cast< const vector< N, T >& >( *( begin() + offset ) );
1266 : }
1267 :
1268 :
1269 : // plane: normal xyz, distance w
1270 : template< size_t M, typename T > template< typename TT >
1271 : inline T vector< M, T >::distance_to_plane( const vector< 3, TT >& point,
1272 : typename enable_if< M == 4, TT >::type* ) const
1273 : {
1274 : const vector< 3, T >& normal = get_sub_vector< 3 >( 0 );
1275 : return normal.dot( point ) + w();
1276 : }
1277 :
1278 :
1279 :
1280 : // plane: normal xyz, distance w
1281 : template< size_t M, typename T >
1282 : template< typename TT >
1283 : vector< 3, T >
1284 : vector< M, T >::project_point_onto_plane( const vector< 3, TT >& point,
1285 : typename enable_if< M == 4, TT >::type* ) const
1286 : {
1287 : const vector< 3, T >& normal = get_sub_vector< 3 >( 0 );
1288 : return point - ( normal * distance_to_plane( point ) );
1289 : }
1290 :
1291 :
1292 :
1293 : template< size_t M, typename T >
1294 23 : bool vector< M, T >::operator==( const vector< M, T >& other ) const
1295 : {
1296 23 : return memcmp( array, other.array, sizeof( array )) == 0;
1297 : }
1298 :
1299 :
1300 : template< size_t M, typename T >
1301 : bool vector< M, T >::operator!=( const vector< M, T >& other ) const
1302 : {
1303 : return ! this->operator==( other );
1304 : }
1305 :
1306 :
1307 : template< size_t M, typename T >
1308 : bool
1309 : vector< M, T >::equals( const vector< M, T >& other, T tolerance ) const
1310 : {
1311 : for( size_t index = 0; index < M; ++index )
1312 : if( fabs( at( index ) - other( index ) ) >= tolerance )
1313 : return false;
1314 : return true;
1315 :
1316 : }
1317 :
1318 :
1319 : template< size_t M, typename T >
1320 : bool
1321 : vector< M, T >::operator<( const vector< M, T >& other ) const
1322 : {
1323 : for(size_t index = 0; index < M; ++index )
1324 : {
1325 : if (at( index ) < other.at( index )) return true;
1326 : if (other.at( index ) < at( index )) return false;
1327 : }
1328 : return false;
1329 : }
1330 :
1331 :
1332 : // to-homogenous-coordinates assignment operator
1333 : // non-chainable because of sfinae
1334 : template< size_t M, typename T > template< size_t N >
1335 : typename enable_if< N == M - 1 >::type*
1336 41 : vector< M, T >::operator=( const vector< N, T >& source_ )
1337 : {
1338 41 : std::copy( source_.begin(), source_.end(), begin() );
1339 41 : at( M - 1 ) = static_cast< T >( 1.0 );
1340 41 : return 0;
1341 : }
1342 :
1343 :
1344 : // from-homogenous-coordinates assignment operator
1345 : // non-chainable because of sfinae
1346 : template< size_t M, typename T > template< size_t N >
1347 : typename enable_if< N == M + 1 >::type*
1348 1 : vector< M, T >::operator=( const vector< N, T >& source_ )
1349 : {
1350 1 : const T w_reci = static_cast< T >( 1.0 ) / source_( M );
1351 1 : iterator it = begin(), it_end = end();
1352 4 : for( size_t index = 0; it != it_end; ++it, ++index )
1353 3 : *it = source_( index ) * w_reci;
1354 1 : return 0;
1355 : }
1356 :
1357 :
1358 : template< size_t M, typename T >
1359 19 : vector< M, T >& vector< M, T >::operator=( const T* c_array )
1360 : {
1361 19 : iter_set( c_array, c_array + M );
1362 19 : return *this;
1363 : }
1364 :
1365 :
1366 :
1367 : template< size_t M, typename T >
1368 : T vector< M, T >::operator=( T filler_value )
1369 : {
1370 : for( size_t index = 0; index < M; ++index )
1371 : at( index ) = filler_value;
1372 : return filler_value;
1373 : }
1374 :
1375 :
1376 :
1377 :
1378 : template< size_t M, typename T >
1379 28 : vector< M, T >& vector< M, T >::operator=( const vector< M, T >& other )
1380 : {
1381 28 : if( this != &other )
1382 28 : memcpy( array, other.array, M * sizeof( T ) );
1383 28 : return *this;
1384 : }
1385 :
1386 :
1387 :
1388 : // returns void to avoid 'silent' loss of precision when chaining
1389 : template< size_t M, typename T > template< typename U >
1390 2 : void vector< M, T >::operator=( const vector< M, U >& source_ )
1391 : {
1392 : typedef typename vector< M, U >::const_iterator u_c_iter;
1393 2 : u_c_iter it = source_.begin(), it_end = source_.end();
1394 8 : for( iterator my_it = begin(); it != it_end; ++it, ++my_it )
1395 6 : *my_it = static_cast< T >( *it );
1396 2 : }
1397 :
1398 :
1399 :
1400 : template< size_t M, typename T >
1401 : template< typename input_iterator_t >
1402 : void
1403 23 : vector< M, T >::iter_set( input_iterator_t begin_, input_iterator_t end_ )
1404 : {
1405 23 : input_iterator_t in_it = begin_;
1406 23 : iterator it = begin(), it_end = end();
1407 110 : for( ; it != it_end && in_it != end_; ++it, ++in_it )
1408 87 : (*it) = static_cast< T >( *in_it );
1409 23 : }
1410 :
1411 : template< size_t M, typename T >
1412 : void vector< M, T >::clamp( const T& min, const T& max )
1413 : {
1414 : for( size_t i = 0; i < M; ++i )
1415 : {
1416 : if( array[i] < min )
1417 : array[i] = min;
1418 : if( array[i] > max )
1419 : array[i] = max;
1420 : }
1421 : }
1422 :
1423 :
1424 :
1425 : template< size_t M, typename T >
1426 : template< typename TT >
1427 : void
1428 : vector< M, T >::scale_to( vector< M, TT >& result_,
1429 : T min_value, T max_value ) const
1430 : {
1431 : T range = max_value-min_value;
1432 : T half_range = range * 0.5;
1433 : T scale = ( 1.0 / range ) * static_cast< T >( std::numeric_limits< TT >::max() );
1434 :
1435 : for( size_t index = 0; index < M; ++index )
1436 : {
1437 : result_.at( index )
1438 : = static_cast< TT >( ( at( index ) + half_range ) * scale );
1439 : }
1440 :
1441 : }
1442 :
1443 :
1444 :
1445 : template< size_t M, typename T >
1446 : inline size_t
1447 : vector< M, T >::size()
1448 : {
1449 : return M;
1450 : }
1451 :
1452 :
1453 :
1454 : template< size_t M, typename T >
1455 : size_t
1456 2 : vector< M, T >::find_min_index() const
1457 : {
1458 2 : return std::min_element( begin(), end() ) - begin();
1459 : }
1460 :
1461 :
1462 :
1463 : template< size_t M, typename T >
1464 : size_t
1465 2 : vector< M, T >::find_max_index() const
1466 : {
1467 2 : return std::max_element( begin(), end() ) - begin();
1468 : }
1469 :
1470 :
1471 :
1472 : template< size_t M, typename T >
1473 : size_t
1474 : vector< M, T >::find_abs_min_index() const
1475 : {
1476 : return std::min_element( begin(), end(), vmml::math::abs_less< T >() ) - begin();
1477 : }
1478 :
1479 :
1480 :
1481 : template< size_t M, typename T >
1482 : size_t
1483 : vector< M, T >::find_abs_max_index() const
1484 : {
1485 : return std::max_element( begin(), end(), vmml::math::abs_greater< T >() ) - begin();
1486 : }
1487 :
1488 :
1489 :
1490 : template< size_t M, typename T >
1491 : T&
1492 2 : vector< M, T >::find_min()
1493 : {
1494 2 : return *std::min_element( begin(), end() );
1495 : }
1496 :
1497 :
1498 :
1499 : template< size_t M, typename T >
1500 : const T&
1501 : vector< M, T >::find_min() const
1502 : {
1503 : return *std::min_element( begin(), end() );
1504 : }
1505 :
1506 :
1507 :
1508 : template< size_t M, typename T >
1509 : T&
1510 2 : vector< M, T >::find_max()
1511 : {
1512 2 : return *std::max_element( begin(), end() );
1513 : }
1514 :
1515 :
1516 :
1517 : template< size_t M, typename T >
1518 : const T&
1519 : vector< M, T >::find_max() const
1520 : {
1521 : return *std::max_element( begin(), end() );
1522 : }
1523 :
1524 :
1525 : template< size_t M, typename T >
1526 : inline typename vector< M, T >::iterator
1527 87 : vector< M, T >::begin()
1528 : {
1529 87 : return array;
1530 : }
1531 :
1532 :
1533 : template< size_t M, typename T >
1534 : inline typename vector< M, T >::iterator
1535 35 : vector< M, T >::end()
1536 : {
1537 35 : return array + M; ;
1538 : }
1539 :
1540 :
1541 : template< size_t M, typename T >
1542 : inline typename vector< M, T >::const_iterator
1543 70 : vector< M, T >::begin() const
1544 : {
1545 70 : return array;
1546 : }
1547 :
1548 :
1549 : template< size_t M, typename T >
1550 : inline typename vector< M, T >::const_iterator
1551 66 : vector< M, T >::end() const
1552 : {
1553 66 : return array + M; ;
1554 : }
1555 :
1556 :
1557 :
1558 : template< size_t M, typename T >
1559 : inline typename vector< M, T >::reverse_iterator
1560 : vector< M, T >::rbegin()
1561 : {
1562 : return array + M - 1;
1563 : }
1564 :
1565 :
1566 : template< size_t M, typename T >
1567 : inline typename vector< M, T >::reverse_iterator
1568 : vector< M, T >::rend()
1569 : {
1570 : return array - 1;
1571 : }
1572 :
1573 :
1574 : template< size_t M, typename T >
1575 : inline typename vector< M, T >::const_reverse_iterator
1576 : vector< M, T >::rbegin() const
1577 : {
1578 : return array + M - 1;
1579 : }
1580 :
1581 :
1582 : template< size_t M, typename T >
1583 : inline typename vector< M, T >::const_reverse_iterator
1584 : vector< M, T >::rend() const
1585 : {
1586 : return array - 1;
1587 : }
1588 :
1589 :
1590 :
1591 : template< size_t M, typename T >
1592 : bool
1593 : vector< M, T >::is_unit_vector() const
1594 : {
1595 : const_iterator it = begin(), it_end = end();
1596 : bool one = false;
1597 : for( ; it != it_end; ++it )
1598 : {
1599 : if ( *it == 1.0 )
1600 : {
1601 : if ( one )
1602 : return false;
1603 : one = true;
1604 : }
1605 : else if ( *it != 0.0 )
1606 : {
1607 : return false;
1608 : }
1609 : }
1610 : return one;
1611 : }
1612 :
1613 :
1614 :
1615 :
1616 : template< size_t M, typename T >
1617 : void
1618 : vector< M, T >::perturb( T perturbation )
1619 : {
1620 : for( iterator it = begin(), it_end = end(); it != it_end; ++it )
1621 : {
1622 : (*it) += ( rand() & 1u ) ? perturbation : -perturbation;
1623 : }
1624 :
1625 : }
1626 :
1627 : template< size_t M, typename T >
1628 : void
1629 1 : vector< M, T >::sqrt_elementwise()
1630 : {
1631 5 : for( iterator it = begin(), it_end = end(); it != it_end; ++it )
1632 : {
1633 4 : (*it) = std::sqrt(*it);
1634 : }
1635 1 : }
1636 :
1637 :
1638 :
1639 : template< size_t M, typename T >
1640 : void
1641 1 : vector< M, T >::reciprocal()
1642 : {
1643 5 : for( iterator it = begin(), it_end = end(); it != it_end; ++it )
1644 : {
1645 4 : (*it) = static_cast< T >( 1.0 ) / (*it);
1646 : }
1647 1 : }
1648 :
1649 :
1650 :
1651 : template< size_t M, typename T >
1652 : void
1653 : vector< M, T >::reciprocal_safe()
1654 : {
1655 : for( iterator it = begin(), it_end = end(); it != it_end; ++it )
1656 : {
1657 : T& v = *it;
1658 :
1659 : if ( v == static_cast< T >( 0 ) )
1660 : v = std::numeric_limits< T >::max();
1661 : else
1662 : v = static_cast< T >( 1.0 ) / v;
1663 : }
1664 : }
1665 :
1666 :
1667 :
1668 : template< size_t M, typename T >
1669 : template< typename TT >
1670 : void
1671 : vector< M, T >::cast_from( const vector< M, TT >& other )
1672 : {
1673 : typedef vmml::vector< M, TT > vector_tt_type ;
1674 : typedef typename vector_tt_type::const_iterator tt_const_iterator;
1675 :
1676 : iterator it = begin(), it_end = end();
1677 : tt_const_iterator other_it = other.begin();
1678 : for( ; it != it_end; ++it, ++other_it )
1679 : {
1680 : *it = static_cast< T >( *other_it );
1681 : }
1682 : }
1683 :
1684 : template< size_t M, typename T >
1685 : size_t
1686 : vector< M, T >::nnz() const
1687 : {
1688 : size_t counter = 0;
1689 :
1690 : const_iterator it = begin(),
1691 : it_end = end();
1692 : for( ; it != it_end; ++it)
1693 : {
1694 : if ( *it != 0 ) {
1695 : ++counter;
1696 : }
1697 : }
1698 :
1699 : return counter;
1700 : }
1701 :
1702 :
1703 : template< size_t M, typename T >
1704 : double
1705 1 : vector< M, T >::norm( ) const
1706 : {
1707 1 : double norm_v = 0.0;
1708 :
1709 1 : const_iterator it = begin(), it_end = end();
1710 5 : for( ; it != it_end; ++it )
1711 : {
1712 4 : norm_v += *it * *it;
1713 : }
1714 :
1715 1 : return std::sqrt(norm_v);
1716 : }
1717 :
1718 : template< size_t M, typename T >
1719 : void
1720 : vector< M, T >::set_random( int seed )
1721 : {
1722 : if ( seed >= 0 )
1723 : srand( seed );
1724 :
1725 : for( size_t i = 0; i < M; ++i )
1726 : {
1727 : const double fillValue = double( rand( )) / double( RAND_MAX );
1728 : at( i ) = -1.0 + 2.0 * fillValue;
1729 : }
1730 : }
1731 :
1732 :
1733 : } // namespace vmml
1734 :
1735 : #endif
|
