vmmlib  1.11.0
Templatized C++ vector and matrix math library
vector.hpp
1 /*
2  * Copyright (c) 2006-2016, 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/enable_if.hpp>
37 
38 #include <algorithm>
39 #include <cmath>
40 #include <cstring>
41 #include <iomanip>
42 #include <iostream>
43 #include <limits>
44 #include <stdexcept>
45 #include <stdint.h>
46 #include <string>
47 #include <vector>
48 
49 namespace vmml
50 {
51 template< size_t M, typename T > class vector
52 {
53 public:
54  typedef T value_type;
55  typedef T* iterator;
56  typedef const T* const_iterator;
57  typedef std::reverse_iterator< iterator > reverse_iterator;
58  typedef std::reverse_iterator< const_iterator > const_reverse_iterator;
59 
60  static const size_t DIMENSION = M;
61 
62  // constructors
63  vector() : array() {} // http://stackoverflow.com/questions/5602030
64  explicit vector( const T& a ); // sets all components to a;
65  vector( const T& x, const T& y );
66  vector( const T& x, const T& y, const T& z );
67  vector( const T& x, const T& y, const T& z, const T& w );
68 
69 #ifndef SWIG
70  // initializes the first M-1 values from vector_, the last from last_
71  vector( const vector< M-1, T >& vector_, T last_ );
72 #endif
73 
74  explicit vector( const T* values );
75 
76 #ifdef __OSG_MATH
77  template< typename OSGVEC3 >
78  explicit vector( const OSGVEC3& from,
79  typename enable_if< M == 3, OSGVEC3 >::type* = 0 );
80 #endif
81 
82  // vec< M > with homogeneous coordinates <-> vec< M-1 > conversion ctor
83  // to-homogenous-coordinates ctor
84  template< size_t N >
85  vector( const vector< N, T >& source_,
86  typename enable_if< N == M - 1 >::type* = 0 );
87 
88  // from-homogenous-coordinates vector
89  template< size_t N >
90  vector( const vector< N, T >& source_,
91  typename enable_if< N == M + 1 >::type* = 0 );
92 
93  template< typename U > vector( const vector< M, U >& source_ );
94 
95  // iterators
96  inline iterator begin();
97  inline iterator end();
98  inline const_iterator begin() const;
99  inline const_iterator end() const;
100  inline reverse_iterator rbegin();
101  inline reverse_iterator rend();
102  inline const_reverse_iterator rbegin() const;
103  inline const_reverse_iterator rend() const;
104 
105 # ifndef VMMLIB_NO_CONVERSION_OPERATORS
106  // conversion operators
107  inline operator T*();
108  inline operator const T*() const;
109 # else
110  inline T& operator[]( size_t index );
111  inline const T& operator[]( size_t index ) const;
112 # endif
113 
114  // accessors
115  inline T& operator()( size_t index );
116  inline const T& operator()( size_t index ) const;
117 
118  inline T& at( size_t index );
119  inline const T& at( size_t index ) const;
120 
121  // element accessors for M <= 4;
122  inline T& x();
123  inline T& y();
124  inline T& z();
125  inline T& w();
126  inline const T& x() const;
127  inline const T& y() const;
128  inline const T& z() const;
129  inline const T& w() const;
130 
131  // pixel color element accessors for M<= 4
132  inline T& r();
133  inline T& g();
134  inline T& b();
135  inline T& a();
136  inline const T& r() const;
137  inline const T& g() const;
138  inline const T& b() const;
139  inline const T& a() const;
140 
141  bool operator==( const vector& other ) const;
142  bool operator!=( const vector& other ) const;
143  bool equals( const vector& other,
144  T tolerance = std::numeric_limits< T >::epsilon( )) const;
145  bool operator<( const vector& other ) const;
146 
147  // remember kids: c_arrays are dangerous and evil!
148  vector& operator=( const T* c_array );
149  T operator=( T filler );
150 
151  vector& operator=( const vector& other );
152 
153  // returns void to avoid 'silent' loss of precision when chaining
154  template< typename U > void operator=( const vector< M, U >& other );
155 
156  // to-homogenous-coordinates assignment operator
157  // non-chainable because of sfinae
158  template< size_t N >
159  typename enable_if< N == M - 1 >::type*
160  operator=( const vector< N, T >& source_ );
161 
162  // from-homogenous-coordinates assignment operator
163  // non-chainable because of sfinae
164  template< size_t N >
165  typename enable_if< N == M + 1 >::type*
166  operator=( const vector< N, T >& source_ );
167 
168  vector operator*( const vector& other ) const;
169  vector operator/( const vector& other ) const;
170  vector operator+( const vector& other ) const;
171  vector operator-( const vector& other ) const;
172 
173  void operator*=( const vector& other );
174  void operator/=( const vector& other );
175  void operator+=( const vector& other );
176  void operator-=( const vector& other );
177 
178  vector operator*( const T other ) const;
179  vector operator/( const T other ) const;
180  vector operator+( const T other ) const;
181  vector operator-( const T other ) const;
182 
183  void operator*=( const T other );
184  void operator/=( const T other );
185  void operator+=( const T other );
186  void operator-=( const T other );
187 
188  vector operator-() const;
189 
190  const vector& negate();
191 
192  void set( T a ); // sets all components to a;
193 #ifndef SWIG
194  void set( const vector< M-1, T >& v, T a );
195 #endif
196  template< size_t N >
197  void set( const vector< N, T >& v );
198 
199  // sets the first few components to a certain value
200  void set( T x, T y );
201  void set( T x, T y, T z );
202  void set( T x, T y, T z, T w );
203 
204  template< typename input_iterator_t >
205  void iter_set( input_iterator_t begin_, input_iterator_t end_ );
206 
207  // compute the cross product of two vectors
208  // note: there's also a free function:
209  // vector<> cross( const vector<>, const vector<> )
210  template< typename TT >
211  vector< M, T >& cross( const vector< M, TT >& b,
212  typename enable_if< M == 3, TT >::type* = 0 );
213 
214  // compute the dot product of two vectors
215  // note: there's also a free function:
216  // T dot( const vector<>, const vector<> );
217  inline T dot( const vector& other ) const;
218 
219  // normalize the vector
220  // note: there's also a free function:
221  // vector<> normalize( const vector<> );
222  inline T normalize();
223 
224  //sets all vector components to random values
225  //remember to set srand( seed );
226  //if seed is set to -1, srand( seed ) was set outside set_random
227  //otherwise srand( seed ) will be called with the given seed
228  void set_random( int seed = -1 );
229 
230  inline T length() const;
231  inline T squared_length() const;
232 
233  inline T distance( const vector& other ) const;
234  inline T squared_distance( const vector& other ) const;
235 
237  T product() const;
238 
239  template< typename TT >
240  vector< 3, T >& rotate( T theta, vector< M, TT > axis,
241  typename enable_if< M == 3, TT >::type* = 0 );
242 
244  template< size_t N, size_t O >
245  vector< N, T > get_sub_vector( typename enable_if< M >= N+O >::type* = 0 )
246  const;
247 
249  template< size_t N, size_t O >
250  void set_sub_vector( const vector< N, T >& sub,
251  typename enable_if< M >= N+O >::type* = 0 );
252 
253  // sphere functions - sphere layout: center xyz, radius w
254  template< typename TT >
255  inline vector< 3, T > project_point_onto_sphere(
256  const vector< 3, TT >& point,
257  typename enable_if< M == 4, TT >::type* = 0 ) const;
258 
259  // returns a negative distance if the point lies in the sphere
260  template< typename TT >
261  inline T distance_to_sphere( const vector< 3, TT >& point,
262  typename enable_if< M == 4, TT >::type* = 0 ) const;
263 
264  // plane functions - plane layout; normal xyz, distance w
265  template< typename TT >
266  inline T distance_to_plane( const vector< 3, TT >& point,
267  typename enable_if< M == 4, TT >::type* = 0 ) const;
268 
269  template< typename TT >
270  inline vector< 3, T > project_point_onto_plane(
271  const vector< 3, TT >& point,
272  typename enable_if< M == 4, TT >::type* = 0 ) const;
273 
274  // returns the index of the minimal resp. maximal value in the vector
275  size_t find_min_index() const;
276  size_t find_max_index() const;
277 
278  // returns minimal resp. maximal value in the vector
279  T& find_min();
280  T& find_max();
281  const T& find_min() const;
282  const T& find_max() const;
283 
284  void clamp( const T& min = 0.0, const T& max = 1.0 );
285 
286  inline static size_t size(); // returns M
287 
288  bool is_unit_vector() const;
289 
290  // perturbs each component by randomly + or - the perturbation parameter
291  void perturb( T perturbation = 0.0001 );
292 
293  void sqrt_elementwise();
294  double norm() const; //l2 norm
295 
296  // computes the reciprocal value for each component, x = 1/x;
297  // WARNING: might result in nans if division by 0!
298  void reciprocal();
299  // computes the reciprocal value for each component, x = 1/x;
300  // checks every component for 0, sets to max value if zero.
301  void reciprocal_safe();
302 
303  template< typename TT >
304  void cast_from( const vector< M, TT >& other );
305 
306  size_t nnz() const;
307 
308  // test each component of the vector for isnan and isinf
309  // inline bool is_valid() const; -> moved to class validator
310 
311  friend std::ostream& operator<< ( std::ostream& os, const vector& vector_ )
312  {
313  const std::ios::fmtflags flags = os.flags();
314  const int prec = os.precision();
315 
316  os.setf( std::ios::right, std::ios::adjustfield );
317  os.precision( 5 );
318  os << "[ ";
319  for( size_t index = 0; index < M; ++index )
320  os << std::setw(10) << vector_.at( index ) << " ";
321  os << "]";
322  os.precision( prec );
323  os.setf( flags );
324  return os;
325  }
326 
327  T array[ M ];
328 
329 #ifndef SWIG
330  // Vector3 defaults
331  static const vector FORWARD;
332  static const vector BACKWARD;
333  static const vector UP;
334  static const vector DOWN;
335  static const vector LEFT;
336  static const vector RIGHT;
337 
338  static const vector ONE;
339  static const vector ZERO;
340 
341  // Unit vectors
342  static const vector UNIT_X;
343  static const vector UNIT_Y;
344  static const vector UNIT_Z;
345 #endif
346 
347 }; // class vector
348 
349 //
350 // typedefs and statics
351 //
352 #ifndef SWIG
353 template< size_t M, typename T >
354 const vector< M, T > vector< M, T >::FORWARD( 0, 0, -1 );
355 template< size_t M, typename T >
357 template< size_t M, typename T >
358 const vector< M, T > vector< M, T >::UP( 0, 1, 0 );
359 template< size_t M, typename T >
360 const vector< M, T > vector< M, T >::DOWN( 0, -1, 0 );
361 template< size_t M, typename T >
362 const vector< M, T > vector< M, T >::LEFT( -1, 0, 0 );
363 template< size_t M, typename T >
364 const vector< M, T > vector< M, T >::RIGHT( 1, 0, 0 );
365 template< size_t M, typename T >
366 
367 const vector< M, T > vector< M, T >::ONE( static_cast< T >( 1 ));
368 template< size_t M, typename T >
369 const vector< M, T > vector< M, T >::ZERO( static_cast< T >( 0 ));
370 template< size_t M, typename T >
371 
372 const vector< M, T > vector< M, T >::UNIT_X( 1, 0, 0 );
373 template< size_t M, typename T >
374 const vector< M, T > vector< M, T >::UNIT_Y( 0, 1, 0 );
375 template< size_t M, typename T >
376 const vector< M, T > vector< M, T >::UNIT_Z( 0, 0, 1 );
377 #endif
378 
379 //
380 // some free functions for convenience
381 //
382 
383 template< size_t M, typename T >
384 bool equals( const vector< M, T >& a, const vector< M, T >& b )
385 {
386  return a.equals( b );
387 }
388 
389 // allows float * vector, not only vector * float
390 template< size_t M, typename T >
391 static vector< M, T > operator* ( T factor, const vector< M, T >& vector_ )
392 {
393  return vector_ * factor;
394 }
395 
396 template< size_t M, typename T >
397 inline T dot( const vector< M, T >& first, const vector< M, T >& second )
398 {
399  return first.dot( second );
400 }
401 
402 template< size_t M, typename T >
403 inline vector< M, T > cross( vector< M, T > a, const vector< M, T >& b )
404 {
405  return a.cross( b );
406 }
407 
408 template< size_t M, typename T >
409 vector< M, T > compute_normal( const vector< M, T >& a, const vector< M, T >& b,
410  const vector< M, T >& c )
411 {
412  // right hand system, CCW triangle
413  const vector< M, T > u = b - a;
414  const vector< M, T > v = c - a;
415  vector< M, T > w = cross( u, v );
416  w.normalize();
417  return w;
418 }
419 
420 template< typename T >
421 vector< 3, T > rotate( vector< 3, T > vec, const T theta,
422  const vector< 3, T >& axis )
423 {
424  return vec.rotate( theta, axis );
425 }
426 
427 
428 template< size_t M, typename T >
429 inline vector< M, T > normalize( vector< M, T > vector_ )
430 {
431  vector_.normalize();
432  return vector_;
433 }
434 
435 template< typename T >
436 inline vector< 4, T > compute_plane( const vector< 3, T >& a,
437  const vector< 3, T >& b,
438  const vector< 3, T >& c )
439 {
440  const vector< 3, T > cb = b - c;
441  const vector< 3, T > ba = a - b;
442 
443  vector< 4, T > plane = vector< 4, T >( cross( cb, ba ));
444  plane.normalize();
445  plane.w() = -plane.x() * a.x() - plane.y() * a.y() - plane.z() * a.z();
446  return plane;
447 }
448 
449 template< size_t M, typename T >
450 vector< M, T >::vector( const T& _a )
451 {
452  for( iterator it = begin(), it_end = end(); it != it_end; ++it )
453  {
454  *it = _a;
455  }
456 }
457 
458 template< size_t M, typename T >
459 vector< M, T >::vector( const T& _x, const T& _y )
460 {
461  array[ 0 ] = _x;
462  array[ 1 ] = _y;
463 }
464 
465 template< size_t M, typename T >
466 vector< M, T >::vector( const T& _x, const T& _y, const T& _z )
467 {
468  array[ 0 ] = _x;
469  array[ 1 ] = _y;
470  array[ 2 ] = _z;
471 }
472 
473 template< size_t M, typename T >
474 vector< M, T >::vector( const T& _x, const T& _y, const T& _z, const T& _w )
475 {
476  array[ 0 ] = _x;
477  array[ 1 ] = _y;
478  array[ 2 ] = _z;
479  array[ 3 ] = _w;
480 }
481 
482 template< size_t M, typename T >
483 vector< M, T >::vector( const T* values )
484 {
485  memcpy( array, values, M * sizeof( T ));
486 }
487 
488 #ifdef __OSG_MATH
489 template< size_t M, typename T >
490 template< typename OSGVEC3 >
491 vector< M, T >::vector( const OSGVEC3& from,
492  typename enable_if< M == 3, OSGVEC3 >::type* )
493 {
494  array[ 0 ] = from.x();
495  array[ 1 ] = from.y();
496  array[ 2 ] = from.z();
497 }
498 #endif
499 
500 #ifndef SWIG
501 template< size_t M, typename T >
502 // initializes the first M-1 values from vector_, the last from last_
503 vector< M, T >::vector( const vector< M-1, T >& vector_, T last_ )
504 {
505  typename vector< M-1, T >::const_iterator
506  it = vector_.begin(), it_end = vector_.end();
507 
508  iterator my_it = begin();
509 
510  for( ; it != it_end; ++it, ++my_it )
511  {
512  (*my_it) = *it;
513  }
514  (*my_it) = last_;
515 }
516 #endif
517 
518 // to-homogenous-coordinates ctor
519 template< size_t M, typename T >
520 template< size_t N >
521 vector< M, T >::vector( const vector< N, T >& source_,
522  typename enable_if< N == M - 1 >::type* )
523 {
524  (*this) = source_;
525 }
526 
527 // from-homogenous-coordinates ctor
528 template< size_t M, typename T >
529 template< size_t N >
530 vector< M, T >::vector( const vector< N, T >& source_,
531  typename enable_if< N == M + 1 >::type* )
532 {
533  (*this) = source_;
534 }
535 
536 template< size_t M, typename T >
537 template< typename U >
538 vector< M, T >::vector( const vector< M, U >& source_ )
539 {
540  (*this) = source_;
541 }
542 
543 template< size_t M, typename T > void vector< M, T >::set( T _a )
544 {
545  for( iterator it = begin(), it_end = end(); it != it_end; ++it )
546  *it = _a;
547 }
548 
549 #ifndef SWIG
550 template< size_t M, typename T >
551 void vector< M, T >::set( const vector< M-1, T >& v, T _a )
552 {
553  memcpy( array, v.array, sizeof( T ) * (M-1) );
554  at( M-1 ) = _a;
555 }
556 #endif
557 
558 template< size_t M, typename T > template< size_t N >
559 void vector< M, T >::set( const vector< N, T >& v )
560 {
561  size_t minimum = M;
562  if (N < M) minimum = N;
563  memcpy( array, v.array, sizeof( T ) * minimum );
564 }
565 
566 template< size_t M, typename T >
567 void vector< M, T >::set( T _x, T _y )
568 {
569  array[ 0 ] = _x;
570  array[ 1 ] = _y;
571 }
572 
573 template< size_t M, typename T >
574 void vector< M, T >::set( T _x, T _y, T _z )
575 {
576  array[ 0 ] = _x;
577  array[ 1 ] = _y;
578  array[ 2 ] = _z;
579 }
580 
581 template< size_t M, typename T >
582 void vector< M, T >::set( T _x, T _y, T _z, T _w )
583 {
584  array[ 0 ] = _x;
585  array[ 1 ] = _y;
586  array[ 2 ] = _z;
587  array[ 3 ] = _w;
588 }
589 
590 template< size_t M, typename T >
591 inline T&
592 vector< M, T >::operator()( size_t index )
593 {
594  return at( index );
595 }
596 
597 template< size_t M, typename T >
598 inline const T&
599 vector< M, T >::operator()( size_t index ) const
600 {
601  return at( index );
602 }
603 
604 template< size_t M, typename T >
605 inline T&
606 vector< M, T >::at( size_t index )
607 {
608  if( index >= M )
609  throw std::runtime_error( "at() - index out of bounds" );
610  return array[ index ];
611 }
612 
613 template< size_t M, typename T >
614 inline const T&
615 vector< M, T >::at( size_t index ) const
616 {
617  if ( index >= M )
618  throw std::runtime_error( "at() - index out of bounds" );
619  return array[ index ];
620 }
621 
622 #ifndef VMMLIB_NO_CONVERSION_OPERATORS
623 
624 template< size_t M, typename T >
625 vector< M, T >::operator T*()
626 {
627  return array;
628 }
629 
630 template< size_t M, typename T >
631 vector< M, T >::operator const T*() const
632 {
633  return array;
634 }
635 #else
636 
637 template< size_t M, typename T >
638 T&
639 vector< M, T >::operator[]( size_t index )
640 {
641  return at( index );
642 }
643 
644 template< size_t M, typename T >
645 const T&
646 vector< M, T >::operator[]( size_t index ) const
647 {
648  return at( index );
649 }
650 
651 #endif
652 
653 #if 0
654 template< size_t M, typename T >
655 inline T&
656 vector< M, T >::operator[]( size_t index )
657 {
658  return at( index );
659 }
660 
661 template< size_t M, typename T >
662 inline const T&
663 vector< M, T >::operator[]( size_t index ) const
664 {
665  return at( index );
666 }
667 #endif
668 
669 template< size_t M, typename T >
670 vector< M, T >
671 vector< M, T >::operator*( const vector< M, T >& other ) const
672 {
673  vector< M, T > result;
674  for( size_t index = 0; index < M; ++index )
675  result.at( index ) = at( index ) * other.at( index );
676  return result;
677 }
678 
679 template< size_t M, typename T >
680 vector< M, T >
681 vector< M, T >::operator/( const vector< M, T >& other ) const
682 {
683  vector< M, T > result;
684  for( size_t index = 0; index < M; ++index )
685  result.at( index ) = at( index ) / other.at( index );
686  return result;
687 }
688 
689 template< size_t M, typename T >
690 vector< M, T >
691 vector< M, T >::operator+( const vector< M, T >& other ) const
692 {
693  vector< M, T > result;
694  for( size_t index = 0; index < M; ++index )
695  result.at( index ) = at( index ) + other.at( index );
696  return result;
697 }
698 
699 template< size_t M, typename T >
700 vector< M, T >
701 vector< M, T >::operator-( const vector< M, T >& other ) const
702 {
703  vector< M, T > result;
704  for( size_t index = 0; index < M; ++index )
705  result.at( index ) = at( index ) - other.at( index );
706  return result;
707 }
708 
709 template< size_t M, typename T >
710 void
711 vector< M, T >::operator*=( const vector< M, T >& other )
712 {
713  for( size_t index = 0; index < M; ++index )
714  at( index ) *= other.at( index );
715 }
716 
717 template< size_t M, typename T >
718 void
719 vector< M, T >::operator/=( const vector< M, T >& other )
720 {
721  for( size_t index = 0; index < M; ++index )
722  at( index ) /= other.at( index );
723 }
724 
725 template< size_t M, typename T >
726 void
727 vector< M, T >::operator+=( const vector< M, T >& other )
728 {
729  for( size_t index = 0; index < M; ++index )
730  at( index ) += other.at( index );
731 }
732 
733 template< size_t M, typename T >
734 void
735 vector< M, T >::operator-=( const vector< M, T >& other )
736 {
737  for( size_t index = 0; index < M; ++index )
738  at( index ) -= other.at( index );
739 }
740 
741 template< size_t M, typename T >
742 vector< M, T >
743 vector< M, T >::operator*( const T other ) const
744 {
745  vector< M, T > result;
746  for( size_t index = 0; index < M; ++index )
747  result.at( index ) = at( index ) * other;
748  return result;
749 }
750 
751 template< size_t M, typename T >
752 vector< M, T >
753 vector< M, T >::operator/( const T other ) const
754 {
755  vector< M, T > result;
756  for( size_t index = 0; index < M; ++index )
757  result.at( index ) = at( index ) / other;
758  return result;
759 }
760 
761 template< size_t M, typename T >
762 vector< M, T >
763 vector< M, T >::operator+( const T other ) const
764 {
765  vector< M, T > result;
766  for( size_t index = 0; index < M; ++index )
767  result.at( index ) = at( index ) + other;
768  return result;
769 }
770 
771 template< size_t M, typename T >
772 vector< M, T >
773 vector< M, T >::operator-( const T other ) const
774 {
775  vector< M, T > result;
776  for( size_t index = 0; index < M; ++index )
777  result.at( index ) = at( index ) - other;
778  return result;
779 }
780 
781 template< size_t M, typename T >
782 void
783 vector< M, T >::operator*=( const T other )
784 {
785  for( size_t index = 0; index < M; ++index )
786  at( index ) *= other;
787 }
788 
789 template< size_t M, typename T >
790 void
791 vector< M, T >::operator/=( const T other )
792 {
793  for( size_t index = 0; index < M; ++index )
794  at( index ) /= other;
795 }
796 
797 template< size_t M, typename T >
798 void
799 vector< M, T >::operator+=( const T other )
800 {
801  for( size_t index = 0; index < M; ++index )
802  at( index ) += other;
803 }
804 
805 template< size_t M, typename T >
806 void
807 vector< M, T >::operator-=( const T other )
808 {
809  for( size_t index = 0; index < M; ++index )
810  at( index ) -= other;
811 }
812 
813 template< size_t M, typename T >
814 vector< M, T >
815 vector< M, T >::operator-() const
816 {
817  vector< M, T > v( *this );
818  return v.negate();
819 }
820 
821 template< size_t M, typename T >
822 const vector< M, T >&
823 vector< M, T >::negate()
824 {
825  for( size_t index = 0; index < M; ++index )
826  array[ index ] = -array[ index ];
827  return *this;
828 }
829 
830 template< size_t M, typename T >
831 inline T&
832 vector< M, T >::x()
833 {
834  return array[ 0 ];
835 }
836 
837 template< size_t M, typename T >
838 inline T&
839 vector< M, T >::y()
840 {
841  return array[ 1 ];
842 }
843 
844 template< size_t M, typename T >
845 inline T&
846 vector< M, T >::z()
847 {
848  return array[ 2 ];
849 }
850 
851 template< size_t M, typename T >
852 inline T&
853 vector< M, T >::w()
854 {
855  return array[ 3 ];
856 }
857 
858 template< size_t M, typename T >
859 inline const T&
860 vector< M, T >::x() const
861 {
862  return array[ 0 ];
863 }
864 
865 template< size_t M, typename T >
866 inline const T&
867 vector< M, T >::y() const
868 {
869  return array[ 1 ];
870 }
871 
872 template< size_t M, typename T >
873 inline const T&
874 vector< M, T >::z() const
875 {
876  return array[ 2 ];
877 }
878 
879 template< size_t M, typename T >
880 inline const T&
881 vector< M, T >::w() const
882 {
883  return array[ 3 ];
884 }
885 
886 template< size_t M, typename T >
887 inline T&
888 vector< M, T >::r()
889 {
890  return array[ 0 ];
891 }
892 
893 template< size_t M, typename T >
894 inline T&
895 vector< M, T >::g()
896 {
897  return array[ 1 ];
898 }
899 
900 template< size_t M, typename T >
901 inline T&
902 vector< M, T >::b()
903 {
904  return array[ 2 ];
905 }
906 
907 template< size_t M, typename T >
908 inline T&
909 vector< M, T >::a()
910 {
911  return array[ 3 ];
912 }
913 
914 template< size_t M, typename T >
915 inline const T&
916 vector< M, T >::r() const
917 {
918  return array[ 0 ];
919 }
920 
921 template< size_t M, typename T >
922 inline const T&
923 vector< M, T >::g() const
924 {
925  return array[ 1 ];
926 }
927 
928 template< size_t M, typename T >
929 inline const T&
930 vector< M, T >::b() const
931 {
932  return array[ 2 ];
933 }
934 
935 template< size_t M, typename T >
936 inline const T&
937 vector< M, T >::a() const
938 {
939  return array[ 3 ];
940 }
941 
942 template< size_t M, typename T > template< typename TT >
943 vector< M, T >& vector< M, T >::cross( const vector< M, TT >& rhs,
944  typename enable_if< M == 3, TT >::type* )
945 {
946  const T x_ = y() * rhs.z() - z() * rhs.y();
947  const T y_ = z() * rhs.x() - x() * rhs.z();
948  const T z_ = x() * rhs.y() - y() * rhs.x();
949  x() = x_;
950  y() = y_;
951  z() = z_;
952  return *this;
953 }
954 
955 template< size_t M, typename T >
956 inline T vector< M, T >::dot( const vector< M, T >& other ) const
957 {
958  T tmp = 0.0;
959  for( size_t index = 0; index < M; ++index )
960  tmp += at( index ) * other.at( index );
961 
962  return tmp;
963 }
964 
965 template< size_t M, typename T > inline T vector< M, T >::normalize()
966 {
967  const T len = length();
968  if ( len <= std::numeric_limits< T >::epsilon( ))
969  return 0;
970 
971  const T tmp = 1.0 / len;
972  (*this) *= tmp;
973  return len;
974 }
975 
976 template< size_t M, typename T >
977 inline T vector< M, T >::length() const
978 {
979  return std::sqrt( squared_length() );
980 }
981 
982 template< size_t M, typename T >
983 inline T vector< M, T >::squared_length() const
984 {
985  T _squared_length = 0.0;
986  for( const_iterator it = begin(), it_end = end(); it != it_end; ++it )
987  _squared_length += (*it) * (*it);
988 
989  return _squared_length;
990 }
991 
992 template< size_t M, typename T >
993 inline T
994 vector< M, T >::distance( const vector< M, T >& other ) const
995 {
996  return std::sqrt( squared_distance( other ) );
997 }
998 
999 template< size_t M, typename T >
1000 inline T vector< M, T >::squared_distance( const vector< M, T >& other ) const
1001 {
1002  vector< M, T > tmp( *this );
1003  tmp -= other;
1004  return tmp.squared_length();
1005 }
1006 
1007 template< size_t M, typename T > inline T vector< M, T >::product() const
1008 {
1009  T result = at( 0 );
1010  for( size_t i = 1; i < M; ++i )
1011  result *= at( i );
1012  return result;
1013 }
1014 
1015 template< size_t M, typename T > template< typename TT >
1017  typename enable_if< M==3, TT >::type* )
1018 {
1019  const T costheta = std::cos( theta );
1020  const T sintheta = std::sin( theta );
1021 
1022  axis.normalize();
1023  return *this = vector< 3, T >(
1024  (costheta + ( 1 - costheta ) * axis.x() * axis.x() ) * x() +
1025  (( 1 - costheta ) * axis.x() * axis.y() - axis.z() * sintheta ) * y() +
1026  (( 1 - costheta ) * axis.x() * axis.z() + axis.y() * sintheta ) * z(),
1027 
1028  (( 1 - costheta ) * axis.x() * axis.y() + axis.z() * sintheta ) * x() +
1029  ( costheta + ( 1 - costheta ) * axis.y() * axis.y() ) * y() +
1030  (( 1 - costheta ) * axis.y() * axis.z() - axis.x() * sintheta ) * z(),
1031 
1032  (( 1 - costheta ) * axis.x() * axis.z() - axis.y() * sintheta ) * x() +
1033  (( 1 - costheta ) * axis.y() * axis.z() + axis.x() * sintheta ) * y() +
1034  ( costheta + ( 1 - costheta ) * axis.z() * axis.z() ) * z( ));
1035 }
1036 
1037 // sphere layout: center xyz, radius w
1038 template< size_t M, typename T > template< typename TT > inline vector< 3, T >
1039 vector< M, T >::project_point_onto_sphere( const vector< 3, TT >& point,
1040  typename enable_if< M == 4, TT >::type* ) const
1041 {
1042  const vector< 3, T >& center_ = get_sub_vector< 3 >();
1043 
1044  vector< 3, T > projected_point( point );
1045  projected_point -= center_;
1046  projected_point.normalize();
1047  projected_point *= w();
1048  return center_ + projected_point;
1049 }
1050 
1051 // sphere layout: center xyz, radius w
1052 template< size_t M, typename T > template< typename TT > inline T
1053 vector< M, T >::distance_to_sphere( const vector< 3, TT >& point,
1054  typename enable_if< M == 4, TT >::type* )
1055  const
1056 {
1057  const vector< 3, T >& center_ = get_sub_vector< 3 >();
1058  return ( point - center_ ).length() - w();
1059 }
1060 
1061 template< size_t M, typename T > template< size_t N, size_t O >
1063  typename enable_if< M >= N+O >::type* ) const
1064 {
1065  return vector< N, T >( array + O );
1066 }
1067 
1068 template< size_t M, typename T > template< size_t N, size_t O >
1070  typename enable_if< M >= N+O >::type* )
1071 {
1072  ::memcpy( array + O, sub.array, N * sizeof( T ));
1073 }
1074 
1075 // plane: normal xyz, distance w
1076 template< size_t M, typename T > template< typename TT >
1077 inline T vector< M, T >::distance_to_plane( const vector< 3, TT >& point,
1078  typename enable_if< M == 4, TT >::type* ) const
1079 {
1080  const vector< 3, T >& normal = get_sub_vector< 3 >();
1081  return normal.dot( point ) + w();
1082 }
1083 
1084 // plane: normal xyz, distance w
1085 template< size_t M, typename T > template< typename TT > vector< 3, T >
1086 vector< M, T >::project_point_onto_plane( const vector< 3, TT >& point,
1087  typename enable_if< M == 4, TT >::type* ) const
1088 {
1089  const vector< 3, T >& normal = get_sub_vector< 3 >();
1090  return point - ( normal * distance_to_plane( point ) );
1091 }
1092 
1093 template< size_t M, typename T >
1094 bool vector< M, T >::operator==( const vector< M, T >& other ) const
1095 {
1096  return memcmp( array, other.array, sizeof( array )) == 0;
1097 }
1098 
1099 template< size_t M, typename T >
1100 bool vector< M, T >::operator!=( const vector< M, T >& other ) const
1101 {
1102  return ! this->operator==( other );
1103 }
1104 
1105 template< size_t M, typename T >
1106 bool vector< M, T >::equals( const vector< M, T >& other, T tolerance ) const
1107 {
1108  for( size_t index = 0; index < M; ++index )
1109  if( fabs( at( index ) - other( index ) ) >= tolerance )
1110  return false;
1111  return true;
1112 
1113 }
1114 
1115 template< size_t M, typename T >
1116 bool
1117 vector< M, T >::operator<( const vector< M, T >& other ) const
1118 {
1119  for(size_t index = 0; index < M; ++index )
1120  {
1121  if (at( index ) < other.at( index )) return true;
1122  if (other.at( index ) < at( index )) return false;
1123  }
1124  return false;
1125 }
1126 
1127 // to-homogenous-coordinates assignment operator
1128 // non-chainable because of sfinae
1129 template< size_t M, typename T > template< size_t N >
1130 typename enable_if< N == M - 1 >::type*
1131 vector< M, T >::operator=( const vector< N, T >& source_ )
1132 {
1133  std::copy( source_.begin(), source_.end(), begin() );
1134  at( M - 1 ) = static_cast< T >( 1.0 );
1135  return 0;
1136 }
1137 
1138 // from-homogenous-coordinates assignment operator
1139 // non-chainable because of sfinae
1140 template< size_t M, typename T > template< size_t N >
1141 typename enable_if< N == M + 1 >::type*
1142 vector< M, T >::operator=( const vector< N, T >& source_ )
1143 {
1144  const T w_reci = static_cast< T >( 1.0 ) / source_( M );
1145  iterator it = begin(), it_end = end();
1146  for( size_t index = 0; it != it_end; ++it, ++index )
1147  *it = source_( index ) * w_reci;
1148  return 0;
1149 }
1150 
1151 template< size_t M, typename T >
1152 vector< M, T >& vector< M, T >::operator=( const T* c_array )
1153 {
1154  iter_set( c_array, c_array + M );
1155  return *this;
1156 }
1157 
1158 template< size_t M, typename T >
1159 T vector< M, T >::operator=( T filler_value )
1160 {
1161  for( size_t index = 0; index < M; ++index )
1162  at( index ) = filler_value;
1163  return filler_value;
1164 }
1165 
1166 template< size_t M, typename T >
1167 vector< M, T >& vector< M, T >::operator=( const vector< M, T >& other )
1168 {
1169  if( this != &other )
1170  memcpy( array, other.array, M * sizeof( T ) );
1171  return *this;
1172 }
1173 
1174 // returns void to avoid 'silent' loss of precision when chaining
1175 template< size_t M, typename T > template< typename U >
1176 void vector< M, T >::operator=( const vector< M, U >& source_ )
1177 {
1178  typedef typename vector< M, U >::const_iterator u_c_iter;
1179  u_c_iter it = source_.begin(), it_end = source_.end();
1180  for( iterator my_it = begin(); it != it_end; ++it, ++my_it )
1181  *my_it = static_cast< T >( *it );
1182 }
1183 
1184 template< size_t M, typename T >
1185 template< typename input_iterator_t >
1186 void
1187 vector< M, T >::iter_set( input_iterator_t begin_, input_iterator_t end_ )
1188 {
1189  input_iterator_t in_it = begin_;
1190  iterator it = begin(), it_end = end();
1191  for( ; it != it_end && in_it != end_; ++it, ++in_it )
1192  (*it) = static_cast< T >( *in_it );
1193 }
1194 
1195 template< size_t M, typename T >
1196 void vector< M, T >::clamp( const T& min, const T& max )
1197 {
1198  for( size_t i = 0; i < M; ++i )
1199  {
1200  if( array[i] < min )
1201  array[i] = min;
1202  if( array[i] > max )
1203  array[i] = max;
1204  }
1205 }
1206 
1207 template< size_t M, typename T >
1208 inline size_t
1209 vector< M, T >::size()
1210 {
1211  return M;
1212 }
1213 
1214 template< size_t M, typename T >
1215 size_t
1216 vector< M, T >::find_min_index() const
1217 {
1218  return std::min_element( begin(), end() ) - begin();
1219 }
1220 
1221 template< size_t M, typename T >
1222 size_t
1223 vector< M, T >::find_max_index() const
1224 {
1225  return std::max_element( begin(), end() ) - begin();
1226 }
1227 
1228 template< size_t M, typename T >
1229 T&
1230 vector< M, T >::find_min()
1231 {
1232  return *std::min_element( begin(), end() );
1233 }
1234 
1235 template< size_t M, typename T >
1236 const T&
1237 vector< M, T >::find_min() const
1238 {
1239  return *std::min_element( begin(), end() );
1240 }
1241 
1242 template< size_t M, typename T >
1243 T&
1244 vector< M, T >::find_max()
1245 {
1246  return *std::max_element( begin(), end() );
1247 }
1248 
1249 template< size_t M, typename T >
1250 const T&
1251 vector< M, T >::find_max() const
1252 {
1253  return *std::max_element( begin(), end() );
1254 }
1255 
1256 template< size_t M, typename T >
1257 inline typename vector< M, T >::iterator
1258 vector< M, T >::begin()
1259 {
1260  return array;
1261 }
1262 
1263 template< size_t M, typename T >
1264 inline typename vector< M, T >::iterator
1265 vector< M, T >::end()
1266 {
1267  return array + M; ;
1268 }
1269 
1270 template< size_t M, typename T >
1271 inline typename vector< M, T >::const_iterator
1272 vector< M, T >::begin() const
1273 {
1274  return array;
1275 }
1276 
1277 template< size_t M, typename T >
1278 inline typename vector< M, T >::const_iterator
1279 vector< M, T >::end() const
1280 {
1281  return array + M; ;
1282 }
1283 
1284 template< size_t M, typename T >
1285 inline typename vector< M, T >::reverse_iterator
1286 vector< M, T >::rbegin()
1287 {
1288  return array + M - 1;
1289 }
1290 
1291 template< size_t M, typename T >
1292 inline typename vector< M, T >::reverse_iterator
1293 vector< M, T >::rend()
1294 {
1295  return array - 1;
1296 }
1297 
1298 template< size_t M, typename T >
1299 inline typename vector< M, T >::const_reverse_iterator
1300 vector< M, T >::rbegin() const
1301 {
1302  return array + M - 1;
1303 }
1304 
1305 template< size_t M, typename T >
1306 inline typename vector< M, T >::const_reverse_iterator
1307 vector< M, T >::rend() const
1308 {
1309  return array - 1;
1310 }
1311 
1312 template< size_t M, typename T >
1313 bool
1314 vector< M, T >::is_unit_vector() const
1315 {
1316  const_iterator it = begin(), it_end = end();
1317  bool one = false;
1318  for( ; it != it_end; ++it )
1319  {
1320  if ( *it == 1.0 )
1321  {
1322  if ( one )
1323  return false;
1324  one = true;
1325  }
1326  else if ( *it != 0.0 )
1327  {
1328  return false;
1329  }
1330  }
1331  return one;
1332 }
1333 
1334 template< size_t M, typename T >
1335 void
1336 vector< M, T >::perturb( T perturbation )
1337 {
1338  for( iterator it = begin(), it_end = end(); it != it_end; ++it )
1339  {
1340  (*it) += ( rand() & 1u ) ? perturbation : -perturbation;
1341  }
1342 
1343 }
1344 
1345 template< size_t M, typename T >
1346 void
1347 vector< M, T >::sqrt_elementwise()
1348 {
1349  for( iterator it = begin(), it_end = end(); it != it_end; ++it )
1350  {
1351  (*it) = std::sqrt(*it);
1352  }
1353 }
1354 
1355 template< size_t M, typename T > void vector< M, T >::reciprocal()
1356 {
1357  for( iterator it = begin(), it_end = end(); it != it_end; ++it )
1358  (*it) = static_cast< T >( 1 ) / (*it);
1359 }
1360 
1361 template< size_t M, typename T > void vector< M, T >::reciprocal_safe()
1362 {
1363  for( iterator it = begin(), it_end = end(); it != it_end; ++it )
1364  {
1365  T& v = *it;
1366 
1367  if( v == 0 )
1368  v = std::numeric_limits< T >::max();
1369  else
1370  v = static_cast< T >( 1 ) / v;
1371  }
1372 }
1373 
1374 template< size_t M, typename T >
1375 template< typename TT >
1376 void
1377 vector< M, T >::cast_from( const vector< M, TT >& other )
1378 {
1379  typedef vmml::vector< M, TT > vector_tt_type ;
1380  typedef typename vector_tt_type::const_iterator tt_const_iterator;
1381 
1382  iterator it = begin(), it_end = end();
1383  tt_const_iterator other_it = other.begin();
1384  for( ; it != it_end; ++it, ++other_it )
1385  {
1386  *it = static_cast< T >( *other_it );
1387  }
1388 }
1389 
1390 template< size_t M, typename T >
1391 size_t
1392 vector< M, T >::nnz() const
1393 {
1394  size_t counter = 0;
1395 
1396  const_iterator it = begin(),
1397  it_end = end();
1398  for( ; it != it_end; ++it)
1399  {
1400  if ( *it != 0 ) {
1401  ++counter;
1402  }
1403  }
1404 
1405  return counter;
1406 }
1407 
1408 template< size_t M, typename T >
1409 double
1410 vector< M, T >::norm( ) const
1411 {
1412  double norm_v = 0.0;
1413 
1414  const_iterator it = begin(), it_end = end();
1415  for( ; it != it_end; ++it )
1416  {
1417  norm_v += *it * *it;
1418  }
1419 
1420  return std::sqrt(norm_v);
1421 }
1422 
1423 template< size_t M, typename T >
1424 void
1425 vector< M, T >::set_random( int seed )
1426 {
1427  if ( seed >= 0 )
1428  srand( seed );
1429 
1430  for( size_t i = 0; i < M; ++i )
1431  {
1432  const double fillValue = double( rand( )) / double( RAND_MAX );
1433  at( i ) = -1.0 + 2.0 * fillValue;
1434  }
1435 }
1436 
1437 } // namespace vmml
1438 
1439 #endif
heavily inspired by boost::enable_if http://www.boost.org, file: boost/utility/enable_if.hpp, Copyright 2003 Jaakko Järvi, Jeremiah Willcock, Andrew Lumsdaine
Definition: aabb.hpp:40
T array[M]
storage
Definition: vector.hpp:327
T product() const
Definition: vector.hpp:1007
void set_sub_vector(const vector< N, T > &sub, typename enable_if< M >=N+O >::type *=0)
Set the sub vector of the given length at the given offset.
Definition: vector.hpp:1069
vector< N, T > get_sub_vector(typename enable_if< M >=N+O >::type *=0) const
Definition: vector.hpp:1062