vmmlib  1.10.0
Templatized C++ vector and matrix math library
matrix.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 #ifndef __VMML__MATRIX__HPP__
33 #define __VMML__MATRIX__HPP__
34 
35 #include <vmmlib/enable_if.hpp>
36 #include <vmmlib/math.hpp>
37 #include <vmmlib/matrix_functors.hpp>
38 #include <vmmlib/types.hpp>
39 
40 #include <iostream>
41 #include <iomanip>
42 #include <vector>
43 #include <cstddef>
44 #include <limits>
45 #include <algorithm>
46 #include <string>
47 #include <cstring>
48 #include <fstream> // file I/O
49 
50 namespace vmml
51 {
52 
53 // matrix of type T with m rows and n columns
54 template< size_t M, size_t N, typename T > class matrix
55 {
56 public:
57  typedef T value_type;
58  typedef T* iterator;
59  typedef const T* const_iterator;
60  typedef std::reverse_iterator< iterator > reverse_iterator;
61  typedef std::reverse_iterator< const_iterator > const_reverse_iterator;
62 
63  static const size_t ROWS = M;
64  static const size_t COLS = N;
65 
70  matrix();
71 
76  matrix( const T* begin, const T* end );
77 
78  template< size_t P, size_t Q, typename U >
79  matrix( const matrix< P, Q, U >& source_ );
80 
85  template< size_t O >
86  matrix( const quaternion< T >& rotation, const vector< O, T >& translation,
87  typename enable_if< M == O+1 && N == O+1 && O == 3 >::type* = 0 );
88 
93  template< size_t O >
94  matrix( const vector< O, T >& eye, const vector< O, T >& lookat,
95  const vector< O, T >& up,
96  typename enable_if< M == O+1 && N == O+1 && O == 3 >::type* = 0 );
97 
98  // accessors
99  inline T& operator()( size_t row_index, size_t col_index );
100  inline const T& operator()( size_t row_index, size_t col_index ) const;
101 
102  inline T& at( size_t row_index, size_t col_index );
103  inline const T& at( size_t row_index, size_t col_index ) const;
104 
105  // element iterators - NOTE: column-major order
106  iterator begin();
107  iterator end();
108  const_iterator begin() const;
109  const_iterator end() const;
110 
111  reverse_iterator rbegin();
112  reverse_iterator rend();
113  const_reverse_iterator rbegin() const;
114  const_reverse_iterator rend() const;
115 
116 #ifndef VMMLIB_NO_CONVERSION_OPERATORS
117  // auto conversion operator
118  operator T*();
119  operator const T*() const;
120 #endif
121 
122  bool operator==( const matrix& other ) const;
123  bool operator!=( const matrix& other ) const;
124 
125  // due to limited precision, two 'identical' matrices might seem different.
126  // this function allows to specify a tolerance when comparing matrices.
127  bool equals( const matrix& other,
128  T tolerance = std::numeric_limits< T >::epsilon( )) const;
129 
130  void multiply_piecewise( const matrix& other );
131 
132  // (this) matrix = left matrix_mxp * right matrix_pxn
133  template< size_t P > void multiply( const matrix< M, P, T >& left,
134  const matrix< P, N, T >& right );
135 
136  // convolution operation (extending borders) of (this) matrix and the given kernel
137  template< size_t U, size_t V >
138  void convolve(const matrix< U, V, T >& kernel);
139 
140  // returned matrix_mxp = (this) matrix * other matrix_nxp;
141  // note: using multiply(...) it avoids a copy of the resulting matrix
142  template< size_t P >
143  matrix< M, P, T > operator*( const matrix< N, P, T >& other ) const;
144 
145  // operator *= is only enabled for square matrices
146  template< size_t O, size_t P, typename TT >
147  typename enable_if< M == N && O == P && M == O, TT >::type*
148  operator*=( const matrix< O, P, TT >& right );
149 
150  inline matrix operator+( const matrix& other ) const;
151  inline matrix operator-( const matrix& other ) const;
152 
153  void operator+=( const matrix& other );
154  void operator-=( const matrix& other );
155 
156  void operator+=( T scalar );
157  void operator-=( T scalar );
158 
159  template< size_t O, size_t P, size_t Q, size_t R >
160  typename enable_if< M == O + Q && N == P + R >::type*
161  direct_sum( const matrix< O, P, T >& m0, const matrix< Q, R, T >& m1 );
162 
163  //
164  // matrix-scalar operations / scaling
165  //
166  matrix operator*( T scalar );
167  void operator*=( T scalar );
168 
169  matrix operator/( T scalar );
170  void operator/=( T scalar );
171 
172  //
173  // matrix-vector operations
174  //
175  // transform column vector by matrix ( vec = matrix * vec )
176  vector< M, T > operator*( const vector< N, T >& other ) const;
177 
178  // transform column vector by matrix ( vec = matrix * vec )
179  // assume homogenous coords, e.g. vec3 = mat4x4 * vec3, with w = 1.0
180  template< size_t O >
181  vector< O, T > operator*( const vector< O, T >& vector_ ) const;
182 
183  inline matrix< M, N, T > operator-() const;
184  matrix< M, N, T > negate() const;
185 
186  // compute tensor product: (this) = vector (X) vector
187  void tensor( const vector< M, T >& u, const vector< N, T >& v );
188  // tensor, for api compatibility with old vmmlib version.
189  // WARNING: for M = N = 4 only.
190  template< size_t uM, size_t vM >
191  typename enable_if< uM == 3 && vM == 3 && M == N && M == 4 >::type*
192  tensor( const vector< uM, T >& u, const vector< vM, T >& v );
193 
194  // row_offset and col_offset define the starting indices for the sub_matrix
195  // the sub_matrix is extracted according to the size of the target matrix, i.e. ( OXP )
196  template< size_t O, size_t P >
197  matrix< O, P, T >
198  get_sub_matrix( size_t row_offset, size_t col_offset,
199  typename enable_if< O <= M && P <= N >::type* = 0 ) const;
200 
201  template< size_t O, size_t P >
202  typename enable_if< O <= M && P <= N >::type*
203  get_sub_matrix( matrix< O, P, T >& result,
204  size_t row_offset = 0, size_t col_offset = 0 ) const;
205 
206  template< size_t O, size_t P >
207  typename enable_if< O <= M && P <= N >::type*
208  set_sub_matrix( const matrix< O, P, T >& sub_matrix,
209  size_t row_offset = 0, size_t col_offset = 0 );
210 
211  // copies a transposed version of *this into transposedMatrix
212  void transpose_to( matrix< N, M, T >& transpose_ ) const;
213 
214  //symmetric covariance matrix of a right matrix multiplication: MxN x NxM = MxM
215  void symmetric_covariance( matrix< M, M, T >& cov_m_ ) const;
216 
217 
218  const matrix& operator=( const matrix< M, N, T >& source_ );
219 
220  // these assignment operators return nothing to avoid silent loss of
221  // precision
222  template< size_t P, size_t Q, typename U >
223  void operator=( const matrix< P, Q, U >& source_ );
224  void operator=( const T old_fashioned_matrix[ M ][ N ] );
225 
226  // WARNING: data_array[] must be at least of size M * N - otherwise CRASH!
227  // WARNING: assumes row_by_row layout - if this is not the case,
228  // use matrix::set( data_array, false )
229  void operator=( const T* data_array );
230  void operator=( const std::vector< T >& data );
231 
232  // sets all elements to fill_value
233  void operator=( T fill_value );
234  void fill( T fill_value );
235 
236  // note: this function copies elements until either the matrix is full or
237  // the iterator equals end_.
238  template< typename input_iterator_t >
239  void set( input_iterator_t begin_, input_iterator_t end_,
240  bool row_major_layout = true );
241 
242  //sets all matrix values with random values
243  //remember to set srand( seed );
244  //if seed is set to -1, srand( seed ) was set outside set_random
245  //otherwise srand( seed ) will be called with the given seed
246  void set_random( int seed = -1 );
247 
248  //sets all matrix values with discrete cosine transform coefficients (receive orthonormal coefficients)
249  void set_dct();
250 
251  void zero();
252 
253  double frobenius_norm() const;
254  double p_norm( double p ) const;
255 
256  template< typename TT >
257  void cast_from( const matrix< M, N, TT >& other );
258 
259  void read_csv_file( const std::string& dir_, const std::string& filename_ );
260  void write_csv_file( const std::string& dir_, const std::string& filename_ ) const;
261  void write_to_raw( const std::string& dir_, const std::string& filename_ ) const;
262  void read_from_raw( const std::string& dir_, const std::string& filename_ ) ;
263 
264  template< typename TT >
265  void quantize_to( matrix< M, N, TT >& quantized_, const T& min_value, const T& max_value ) const;
266  template< typename TT >
267  void quantize( matrix< M, N, TT >& quantized_, T& min_value, T& max_value ) const;
268  template< typename TT >
269  void dequantize( matrix< M, N, TT >& quantized_, const TT& min_value, const TT& max_value ) const;
270 
271  void columnwise_sum( vector< N, T>& summed_columns_ ) const;
272  double sum_elements() const;
273 
274  void sum_rows( matrix< M/2, N, T>& other ) const;
275  void sum_columns( matrix< M, N/2, T>& other ) const;
276 
277  template< size_t R >
278  typename enable_if< R == M && R == N >::type*
279  diag( const vector< R, T >& diag_values_ );
280 
281 
282  //Khatri-Rao Product: columns must be of same size
283  template< size_t O >
284  void khatri_rao_product( const matrix< O, N, T >& right_, matrix< M*O, N, T >& prod_ ) const;
285  //Kronecker Product: MxN x_kronecker OxP = M*OxN*P
286  template< size_t O, size_t P >
287  void kronecker_product( const matrix< O, P, T >& right_, matrix< M*O, N*P, T >& result_) const;
288 
289  T get_min() const;
290  T get_max() const;
291  T get_abs_min() const;
292  T get_abs_max() const;
293 
294  //returns number of non-zeros
295  size_t nnz() const;
296  size_t nnz( const T& threshold_ ) const;
297  void threshold( const T& threshold_value_ );
298 
299  vector< M, T > get_column( size_t column_index ) const;
300  void get_column( size_t column_index, vector< M, T>& column ) const;
301 
302  void set_column( size_t index, const vector< M, T >& column );
303 
304  void get_column( size_t index, matrix< M, 1, T >& column ) const;
305  void set_column( size_t index, const matrix< M, 1, T >& column );
306 
307  vector< N, T > get_row( size_t index ) const;
308  void get_row( size_t index, vector< N, T >& row ) const;
309  void set_row( size_t index, const vector< N, T >& row );
310 
311  void get_row( size_t index, matrix< 1, N, T >& row ) const;
312  void set_row( size_t index, const matrix< 1, N, T >& row );
313 
314  size_t size() const; // return M * N;
315 
316  size_t get_number_of_rows() const;
317  size_t get_number_of_columns() const;
318 
319  T det() const;
320 
321  // the return value indicates if the matrix is invertible.
322  // we need a tolerance term since the computation of the determinant is
323  // subject to precision errors.
324  template< size_t O, size_t P, typename TT >
325  bool inverse( matrix< O, P, TT >& inverse_,
326  T tolerance = std::numeric_limits<T>::epsilon(),
327  typename enable_if< M == N && O == P && O == M && M >= 2 && M <= 4, TT >::type* = 0 ) const;
328 
329  template< size_t O, size_t P >
330  typename enable_if< O == P && M == N && O == M && M >= 2 >::type*
331  get_adjugate( matrix< O, P, T >& adjugate ) const;
332 
333  template< size_t O, size_t P >
334  typename enable_if< O == P && M == N && O == M && M >= 2 >::type*
335  get_cofactors( matrix< O, P, T >& cofactors ) const;
336 
337 
338  // returns the determinant of a square matrix M-1, N-1
339  template< size_t O, size_t P >
340  T get_minor( matrix< O, P, T >& minor_, size_t row_to_cut, size_t col_to_cut,
341  typename enable_if< O == M-1 && P == N-1 && M == N && M >= 2 >::type* = 0 ) const;
342 
343  //
344  // 4*4 matrices only
345  //
350  template< typename TT >
351  matrix< M, N, T >& rotate( const TT angle, const vector< M-1, T >& axis,
352  typename enable_if< M == N && M == 4, TT >::type* = 0 );
353 
354  template< typename TT >
355  matrix< M, N, T >& rotate_x( const TT angle,
356  typename enable_if< M == N && M == 4, TT >::type* = 0 );
357 
358  template< typename TT >
359  matrix< M, N, T >& rotate_y( const TT angle,
360  typename enable_if< M == N && M == 4, TT >::type* = 0 );
361 
362  template< typename TT >
363  matrix< M, N, T >& rotate_z( const TT angle,
364  typename enable_if< M == N && M == 4, TT >::type* = 0 );
365 
366  template< typename TT >
367  matrix< M, N, T >& pre_rotate_x( const TT angle,
368  typename enable_if< M == N && M == 4, TT >::type* = 0 );
369 
370  template< typename TT >
371  matrix< M, N, T >& pre_rotate_y( const TT angle,
372  typename enable_if< M == N && M == 4, TT >::type* = 0 );
373 
374  template< typename TT >
375  matrix< M, N, T >& pre_rotate_z( const TT angle,
376  typename enable_if< M == N && M == 4, TT >::type* = 0 );
377 
378  template< typename TT >
379  matrix< M, N, T >& scale( const TT scale[3],
380  typename enable_if< M == N && M == 4, TT >::type* = 0 );
381 
382  template< typename TT >
383  matrix< M, N, T >& scale( const TT x_, const T y, const T z,
384  typename enable_if< M == N && M == 4, TT >::type* = 0 );
385 
386  template< typename TT >
387  matrix< M, N, T >& scale( const vector< 3, TT >& scale_,
388  typename enable_if< M == N && M == 4, TT >::type* = 0 );
389 
390  template< typename TT >
391  matrix< M, N, T >& scale_translation( const TT scale_[3],
392  typename enable_if< M == N && M == 4, TT >::type* = 0 );
393 
394  template< typename TT >
395  matrix< M, N, T >& scale_translation( const vector< 3, TT >& scale_,
396  typename enable_if< M == N && M == 4, TT >::type* = 0 );
397 
398  template< typename TT >
399  matrix< M, N, T >& set_translation( const TT x_, const TT y_, const TT z_,
400  typename enable_if< M == N && M == 4, TT >::type* = 0 );
401 
402  template< typename TT >
403  matrix< M, N, T >& set_translation( const TT trans[3],
404  typename enable_if< M == N && M == 4, TT >::type* = 0 );
405 
406  template< typename TT >
407  matrix< M, N, T >& set_translation( const vector< 3, TT >& t,
408  typename enable_if< M == N && M == 4, TT >::type* = 0 );
409 
410  template< typename TT >
411  void get_translation( vector< N-1, TT >& translation_ ) const;
412 
413  vector< N-1, T > get_translation() const;
414 
419  template< size_t O >
420  void getLookAt( vector< O, T >& eye, vector< O, T >& lookAt,
421  vector< O, T >& up,
422  typename enable_if< M == O+1 && N == O+1 && O == 3 >::type* = 0 ) const;
423 
424  // hack for static-member-init
425  template< typename init_functor_t >
426  static const matrix get_initialized_matrix();
427 
428  inline T& x();
429  inline T& y();
430  inline T& z();
431 
432 
433  // tests every component for isnan && isinf
434  // inline bool is_valid() const; -> moved to class validator
435 
436  // legacy/compatibility accessor
437  struct row_accessor
438  {
439  row_accessor( T* array_ ) : array( array_ ) {}
440  T& operator[]( const size_t col_index )
441  {
442  if ( col_index >= N )
443  throw std::runtime_error( "column index out of bounds" );
444  return array[ col_index * M ];
445  }
446 
447  const T& operator[]( const size_t col_index ) const
448  {
449  if ( col_index >= N )
450  throw std::runtime_error( "column index out of bounds" );
451  return array[ col_index * M ];
452  }
453 
454  private:
455  row_accessor() {} // disallow std ctor
456  T* array;
457  };
458  // this is a hack to allow array-style access to matrix elements
459  // usage: matrix< 2, 2, float > m; m[ 1 ][ 0 ] = 37.0f;
460  inline row_accessor operator[]( const size_t row_index )
461  {
462  if ( row_index > M )
463  throw std::runtime_error( "row index out of bounds" );
464  return row_accessor( array + row_index );
465  }
466 
467  // this is a hack to remove a warning about implicit conversions
468  inline row_accessor operator[]( int row_index )
469  {
470  return ( *this )[ size_t ( row_index ) ];
471  }
472 
473  friend std::ostream& operator << ( std::ostream& os,
474  const matrix< M, N, T >& matrix )
475  {
476  const std::ios::fmtflags flags = os.flags();
477  const int prec = os.precision();
478 
479  os.setf( std::ios::right, std::ios::adjustfield );
480  os.precision( 5 );
481 
482  for( size_t row_index = 0; row_index < M; ++row_index )
483  {
484  os << "|";
485  for( size_t col_index = 0; col_index < N; ++col_index )
486  {
487  os << std::setw(10) << matrix.at( row_index, col_index ) << " ";
488  }
489  os << "|" << std::endl;
490  }
491  os.precision( prec );
492  os.setf( flags );
493  return os;
494  };
495 
496  // protected:
497  T array[ M * N ];
498 
499  // static members
500  static const matrix< M, N, T > IDENTITY;
501  static const matrix< M, N, T > ZERO;
502 
503 }; // class matrix
504 }
505 
506 #include <vmmlib/quaternion.hpp>
507 #include <vmmlib/vector.hpp>
508 
509 namespace vmml
510 {
511 /*
512  * free functions
513  */
514 template< size_t M, size_t N, typename T >
515 inline void multiply( const matrix< M, N, T >& left,
516  const matrix< M, N, T >& right,
517  matrix< M, N, T >& result )
518 {
519  result.multiply( left, right );
520 }
521 
522 
523 
524 template< size_t M, size_t N, typename T >
525 template< size_t U, size_t V >
526 void matrix< M, N, T>::convolve(const matrix< U, V, T >& kernel)
527 {
528  matrix< M, N, T> temp; // do not override original values instantly as old values are needed for calculation
529 
530  for(size_t y_ = 0; y_ < N; ++y_)
531  {
532  for(size_t x_ = 0; x_ < M; ++x_)
533  {
534  double sum = 0.0;
535 
536  for(size_t j = 0; j < V; ++j)
537  {
538  int srcy = y_ - V/2 + j;
539 
540  // Extending border values
541  if(srcy < 0) srcy = 0;
542  if(srcy >= int(N)) srcy = N-1;
543 
544  for(size_t i = 0; i < U; ++i)
545  {
546  int srcx = x_ - U/2 + i;
547 
548  // Extending border values
549  if(srcx < 0) srcx = 0;
550  if(srcx >= int(M)) srcx = M-1;
551 
552  sum += kernel.at(j,i) * at(srcy,srcx);
553  }
554  }
555  temp.at(y_,x_) = sum;
556  }
557  }
558 
559  *this = temp;
560 }
561 
562 
563 
564 template< size_t M, size_t N, size_t P, typename T >
565 inline void
566 multiply( const matrix< M, P, T >& left, const matrix< P, N, T >& right,
567  matrix< M, N, T >& result )
568 {
569  result.multiply( left, right );
570 }
571 
572 
573 template< size_t M, size_t N, typename T >
574 inline typename enable_if< M == N >::type* identity( matrix< M, N, T >& m )
575 {
576  m = static_cast< T >( 0.0 );
577  for( size_t index = 0; index < M; ++index )
578  m( index, index ) = static_cast< T >( 1.0 );
579  return 0; // for sfinae
580 }
581 
582 
583 
584 template< typename T >
585 inline T compute_determinant( const matrix< 1, 1, T >& matrix_ )
586 {
587  return matrix_.array[ 0 ];
588 }
589 
590 
591 
592 template< typename T >
593 inline T compute_determinant( const matrix< 2, 2, T >& matrix_ )
594 {
595  const T& a = matrix_( 0, 0 );
596  const T& b = matrix_( 0, 1 );
597  const T& c = matrix_( 1, 0 );
598  const T& d = matrix_( 1, 1 );
599  return a * d - b * c;
600 }
601 
602 
603 
604 template< typename T >
605 inline T compute_determinant( const matrix< 3, 3, T >& m_ )
606 {
607  return
608  m_( 0,0 ) * ( m_( 1,1 ) * m_( 2,2 ) - m_( 1,2 ) * m_( 2,1 ) )
609  + m_( 0,1 ) * ( m_( 1,2 ) * m_( 2,0 ) - m_( 1,0 ) * m_( 2,2 ) )
610  + m_( 0,2 ) * ( m_( 1,0 ) * m_( 2,1 ) - m_( 1,1 ) * m_( 2,0 ) );
611 }
612 
613 
614 template< typename T >
615 inline T compute_determinant( const matrix< 4, 4, T >& m )
616 {
617  T m00 = m( 0, 0 );
618  T m10 = m( 1, 0 );
619  T m20 = m( 2, 0 );
620  T m30 = m( 3, 0 );
621  T m01 = m( 0, 1 );
622  T m11 = m( 1, 1 );
623  T m21 = m( 2, 1 );
624  T m31 = m( 3, 1 );
625  T m02 = m( 0, 2 );
626  T m12 = m( 1, 2 );
627  T m22 = m( 2, 2 );
628  T m32 = m( 3, 2 );
629  T m03 = m( 0, 3 );
630  T m13 = m( 1, 3 );
631  T m23 = m( 2, 3 );
632  T m33 = m( 3, 3 );
633 
634  return
635  m03 * m12 * m21 * m30
636  - m02 * m13 * m21 * m30
637  - m03 * m11 * m22 * m30
638  + m01 * m13 * m22 * m30
639  + m02 * m11 * m23 * m30
640  - m01 * m12 * m23 * m30
641  - m03 * m12 * m20 * m31
642  + m02 * m13 * m20 * m31
643  + m03 * m10 * m22 * m31
644  - m00 * m13 * m22 * m31
645  - m02 * m10 * m23 * m31
646  + m00 * m12 * m23 * m31
647  + m03 * m11 * m20 * m32
648  - m01 * m13 * m20 * m32
649  - m03 * m10 * m21 * m32
650  + m00 * m13 * m21 * m32
651  + m01 * m10 * m23 * m32
652  - m00 * m11 * m23 * m32
653  - m02 * m11 * m20 * m33
654  + m01 * m12 * m20 * m33
655  + m02 * m10 * m21 * m33
656  - m00 * m12 * m21 * m33
657  - m01 * m10 * m22 * m33
658  + m00 * m11 * m22 * m33;
659 }
660 
661 
662 
663 template< typename T >
664 bool compute_inverse( const matrix< 2, 2, T >& m_, matrix< 2, 2, T >& inverse_,
665  T tolerance_ = std::numeric_limits<T>::epsilon())
666 {
667  const T det_ = compute_determinant( m_ );
668  if( fabs( det_ ) < tolerance_ )
669  return false;
670 
671  const T detinv = static_cast< T >( 1.0 ) / det_;
672 
673  m_.get_adjugate( inverse_ ); // set inverse_ to the adjugate of M
674  inverse_ *= detinv;
675  return true;
676 }
677 
678 
679 
680 template< typename T >
681 bool compute_inverse( const matrix< 3, 3, T >& m_, matrix< 3, 3, T >& inverse_,
682  T tolerance_ = std::numeric_limits<T>::epsilon() )
683 {
684  // Invert a 3x3 using cofactors. This is about 8 times faster than
685  // the Numerical Recipes code which uses Gaussian elimination.
686  inverse_.at( 0, 0 ) = m_.at( 1, 1 ) * m_.at( 2, 2 ) - m_.at( 1, 2 ) * m_.at( 2, 1 );
687  inverse_.at( 0, 1 ) = m_.at( 0, 2 ) * m_.at( 2, 1 ) - m_.at( 0, 1 ) * m_.at( 2, 2 );
688  inverse_.at( 0, 2 ) = m_.at( 0, 1 ) * m_.at( 1, 2 ) - m_.at( 0, 2 ) * m_.at( 1, 1 );
689  inverse_.at( 1, 0 ) = m_.at( 1, 2 ) * m_.at( 2, 0 ) - m_.at( 1, 0 ) * m_.at( 2, 2 );
690  inverse_.at( 1, 1 ) = m_.at( 0, 0 ) * m_.at( 2, 2 ) - m_.at( 0, 2 ) * m_.at( 2, 0 );
691  inverse_.at( 1, 2 ) = m_.at( 0, 2 ) * m_.at( 1, 0 ) - m_.at( 0, 0 ) * m_.at( 1, 2 );
692  inverse_.at( 2, 0 ) = m_.at( 1, 0 ) * m_.at( 2, 1 ) - m_.at( 1, 1 ) * m_.at( 2, 0 );
693  inverse_.at( 2, 1 ) = m_.at( 0, 1 ) * m_.at( 2, 0 ) - m_.at( 0, 0 ) * m_.at( 2, 1 );
694  inverse_.at( 2, 2 ) = m_.at( 0, 0 ) * m_.at( 1, 1 ) - m_.at( 0, 1 ) * m_.at( 1, 0 );
695 
696  const T determinant =
697  m_.at( 0, 0 ) * inverse_.at( 0, 0 )
698  + m_.at( 0, 1 ) * inverse_.at( 1, 0 )
699  + m_.at( 0, 2 ) * inverse_.at( 2, 0 );
700 
701  if ( fabs( determinant ) <= tolerance_ )
702  return false; // matrix is not invertible
703 
704  const T detinv = static_cast< T >( 1.0 ) / determinant;
705 
706  inverse_.at( 0, 0 ) *= detinv;
707  inverse_.at( 0, 1 ) *= detinv;
708  inverse_.at( 0, 2 ) *= detinv;
709  inverse_.at( 1, 0 ) *= detinv;
710  inverse_.at( 1, 1 ) *= detinv;
711  inverse_.at( 1, 2 ) *= detinv;
712  inverse_.at( 2, 0 ) *= detinv;
713  inverse_.at( 2, 1 ) *= detinv;
714  inverse_.at( 2, 2 ) *= detinv;
715 
716  return true;
717 }
718 
719 
720 template< typename T >
721 bool compute_inverse( const matrix< 4, 4, T >& m_,
722  matrix< 4, 4, T >& inv_,
723  T tolerance_ = std::numeric_limits<T>::epsilon() )
724 {
725  const T* array = m_.array;
726 
727  // tuned version from Claude Knaus
728  /* first set of 2x2 determinants: 12 multiplications, 6 additions */
729  const T t1[6] = { array[ 2] * array[ 7] - array[ 6] * array[ 3],
730  array[ 2] * array[11] - array[10] * array[ 3],
731  array[ 2] * array[15] - array[14] * array[ 3],
732  array[ 6] * array[11] - array[10] * array[ 7],
733  array[ 6] * array[15] - array[14] * array[ 7],
734  array[10] * array[15] - array[14] * array[11] };
735 
736  /* first half of comatrix: 24 multiplications, 16 additions */
737  inv_.array[0] = array[ 5] * t1[5] - array[ 9] * t1[4] + array[13] * t1[3];
738  inv_.array[1] = array[ 9] * t1[2] - array[13] * t1[1] - array[ 1] * t1[5];
739  inv_.array[2] = array[13] * t1[0] - array[ 5] * t1[2] + array[ 1] * t1[4];
740  inv_.array[3] = array[ 5] * t1[1] - array[ 1] * t1[3] - array[ 9] * t1[0];
741  inv_.array[4] = array[ 8] * t1[4] - array[ 4] * t1[5] - array[12] * t1[3];
742  inv_.array[5] = array[ 0] * t1[5] - array[ 8] * t1[2] + array[12] * t1[1];
743  inv_.array[6] = array[ 4] * t1[2] - array[12] * t1[0] - array[ 0] * t1[4];
744  inv_.array[7] = array[ 0] * t1[3] - array[ 4] * t1[1] + array[ 8] * t1[0];
745 
746  /* second set of 2x2 determinants: 12 multiplications, 6 additions */
747  const T t2[6] = { array[ 0] * array[ 5] - array[ 4] * array[ 1],
748  array[ 0] * array[ 9] - array[ 8] * array[ 1],
749  array[ 0] * array[13] - array[12] * array[ 1],
750  array[ 4] * array[ 9] - array[ 8] * array[ 5],
751  array[ 4] * array[13] - array[12] * array[ 5],
752  array[ 8] * array[13] - array[12] * array[ 9] };
753 
754  /* second half of comatrix: 24 multiplications, 16 additions */
755  inv_.array[8] = array[ 7] * t2[5] - array[11] * t2[4] + array[15] * t2[3];
756  inv_.array[9] = array[11] * t2[2] - array[15] * t2[1] - array[ 3] * t2[5];
757  inv_.array[10] = array[15] * t2[0] - array[ 7] * t2[2] + array[ 3] * t2[4];
758  inv_.array[11] = array[ 7] * t2[1] - array[ 3] * t2[3] - array[11] * t2[0];
759  inv_.array[12] = array[10] * t2[4] - array[ 6] * t2[5] - array[14] * t2[3];
760  inv_.array[13] = array[ 2] * t2[5] - array[10] * t2[2] + array[14] * t2[1];
761  inv_.array[14] = array[ 6] * t2[2] - array[14] * t2[0] - array[ 2] * t2[4];
762  inv_.array[15] = array[ 2] * t2[3] - array[ 6] * t2[1] + array[10] * t2[0];
763 
764  /* determinant: 4 multiplications, 3 additions */
765  const T determinant = array[0] * inv_.array[0] + array[4] * inv_.array[1] +
766  array[8] * inv_.array[2] + array[12] * inv_.array[3];
767 
768  if( fabs( determinant ) <= tolerance_ )
769  return false; // matrix is not invertible
770 
771  /* division: 16 multiplications, 1 division */
772  const T detinv = static_cast< T >( 1.0 ) / determinant;
773 
774  for( size_t i = 0; i != 16; ++i )
775  inv_.array[i] *= detinv;
776 
777  return true;
778 }
779 
780 
781 // this function returns the transpose of a matrix
782 // however, using matrix::transpose_to( .. ) avoids the copy.
783 template< size_t M, size_t N, typename T >
784 inline matrix< N, M, T >
785 transpose( const matrix< M, N, T >& matrix_ )
786 {
787  matrix< N, M, T > transpose_;
788  matrix_.transpose_to( transpose_ );
789  return transpose_;
790 }
791 
792 
793 template< size_t M, size_t N, typename T >
794 bool is_positive_definite( const matrix< M, N, T >& matrix_,
795  const T limit = -1e12,
796  typename enable_if< M == N && M <= 3 >::type* = 0 )
797 {
798  if ( matrix_.at( 0, 0 ) < limit )
799  return false;
800 
801  // sylvester criterion
802  if ( M > 1 )
803  {
804  matrix< 2, 2, T > m;
805  matrix_.get_sub_matrix( m, 0, 0 );
806  if ( compute_determinant( m ) < limit )
807  return false;
808  }
809 
810  if ( M > 2 )
811  {
812  matrix< 3, 3, T > m;
813  matrix_.get_sub_matrix( m, 0, 0 );
814  if ( compute_determinant( m ) < limit )
815  return false;
816  }
817 
818  return true;
819 }
820 
821 
822 template< size_t M, size_t N, typename T >
823 matrix< M, N, T >::matrix()
824  : array() // http://stackoverflow.com/questions/5602030
825 {
826  if( M == N )
827  for( size_t i = 0; i < M; ++i )
828  at( i, i ) = static_cast< T >( 1.0 );
829 }
830 
831 template< size_t M, size_t N, typename T >
832 matrix< M, N, T >::matrix( const T* begin_, const T* end_ )
833  : array() // http://stackoverflow.com/questions/5602030
834 {
835  const T* to = std::min( end_, begin_ + M*N );
836  ::memcpy( array, begin_, (to - begin_) * sizeof( T ));
837 }
838 
839 template< size_t M, size_t N, typename T >
840 template< size_t P, size_t Q, typename U >
841 matrix< M, N, T >::matrix( const matrix< P, Q, U >& source_ )
842 {
843  (*this) = source_;
844 }
845 
846 template< size_t M, size_t N, typename T > template< size_t O >
847 matrix< M, N, T >::matrix( const quaternion< T >& rotation,
848  const vector< O, T >& translation,
849  typename enable_if< M == O+1 && N == O+1 && O == 3 >::type* )
850 {
851  rotation.get_rotation_matrix( *this );
852  set_translation( translation );
853  at( 3, 0 ) = 0;
854  at( 3, 1 ) = 0;
855  at( 3, 2 ) = 0;
856  at( 3, 3 ) = 1;
857 }
858 
859 template< size_t M, size_t N, typename T >
860 template< size_t O >
861 matrix< M, N, T >::matrix( const vector< O, T >& eye,
862  const vector< O, T >& lookat,
863  const vector< O, T >& up,
864  typename enable_if< M == O+1 && N == O+1 && O == 3 >::type* )
865 {
866  *this = matrix< 4, 4, T >::IDENTITY;
867 
868  const vector< 3, T > f( vmml::normalize( lookat - eye ));
869  const vector< 3, T > s( vmml::normalize( vmml::cross( f, up )));
870  const vector< 3, T > u( vmml::cross( s, f ));
871 
872  at( 0, 0 ) = s.x();
873  at( 0, 1 ) = s.y();
874  at( 0, 2 ) = s.z();
875  at( 1, 0 ) = u.x();
876  at( 1, 1 ) = u.y();
877  at( 1, 2 ) = u.z();
878  at( 2, 0 ) = -f.x();
879  at( 2, 1 ) = -f.y();
880  at( 2, 2 ) = -f.z();
881  at( 0, 3 ) = -vmml::dot( s, eye );
882  at( 1, 3 ) = -vmml::dot( u, eye );
883  at( 2, 3 ) = vmml::dot( f, eye );
884 }
885 
886 template< size_t M, size_t N, typename T >
887 inline T& matrix< M, N, T >::at( size_t row_index, size_t col_index )
888 {
889 #ifdef VMMLIB_SAFE_ACCESSORS
890  if ( row_index >= M || col_index >= N )
891  throw std::runtime_error( "at( row, col ) - index out of bounds" );
892 #endif
893  return array[ col_index * M + row_index ];
894 }
895 
896 
897 
898 template< size_t M, size_t N, typename T >
899 const inline T&
900 matrix< M, N, T >::at( size_t row_index, size_t col_index ) const
901 {
902 #ifdef VMMLIB_SAFE_ACCESSORS
903  if ( row_index >= M || col_index >= N )
904  throw std::runtime_error( "at( row, col ) - index out of bounds" );
905 #endif
906  return array[ col_index * M + row_index ];
907 }
908 
909 
910 template< size_t M, size_t N, typename T >
911 inline T&
912 matrix< M, N, T >::operator()( size_t row_index, size_t col_index )
913 {
914  return at( row_index, col_index );
915 }
916 
917 
918 
919 template< size_t M, size_t N, typename T >
920 const inline T&
921 matrix< M, N, T >::operator()( size_t row_index, size_t col_index ) const
922 {
923  return at( row_index, col_index );
924 }
925 
926 #ifndef VMMLIB_NO_CONVERSION_OPERATORS
927 
928 template< size_t M, size_t N, typename T >
929 matrix< M, N, T >::
930 operator T*()
931 {
932  return array;
933 }
934 
935 
936 
937 template< size_t M, size_t N, typename T >
938 matrix< M, N, T >::
939 operator const T*() const
940 {
941  return array;
942 }
943 
944 #endif
945 
946 
947 template< size_t M, size_t N, typename T >
948 bool
949 matrix< M, N, T >::
950 operator==( const matrix< M, N, T >& other ) const
951 {
952  bool is_ok = true;
953  for( size_t i = 0; is_ok && i < M * N; ++i )
954  {
955  is_ok = array[ i ] == other.array[ i ];
956  }
957  return is_ok;
958 }
959 
960 
961 
962 
963 template< size_t M, size_t N, typename T >
964 bool
965 matrix< M, N, T >::
966 operator!=( const matrix< M, N, T >& other ) const
967 {
968  return ! operator==( other );
969 }
970 
971 
972 
973 template< size_t M, size_t N, typename T >
974 bool matrix< M, N, T >::equals( const matrix< M, N, T >& other,
975  const T tolerance ) const
976 {
977  for( size_t row_index = 0; row_index < M; row_index++)
978  for( size_t col_index = 0; col_index < N; col_index++)
979  if( std::abs( at( row_index, col_index ) -
980  other( row_index, col_index )) > tolerance )
981  {
982  return false;
983  }
984  return true;
985 }
986 
987 template< size_t M, size_t N, typename T >
988 const matrix< M, N, T >&
989 matrix< M, N, T >::operator=( const matrix< M, N, T >& source_ )
990 {
991  memcpy( array, source_.array, M * N * sizeof( T ));
992  return *this;
993 }
994 
995 template< size_t M, size_t N, typename T >
996 template< size_t P, size_t Q, typename U >
997 void matrix< M, N, T >::operator=( const matrix< P, Q, U >& source_ )
998 {
999  const size_t minL = P < M ? P : M;
1000  const size_t minC = Q < N ? Q : N;
1001 
1002  if( minL < M || minC < N )
1003  zero();
1004 
1005  for ( size_t i = 0 ; i < minL ; i++ )
1006  for ( size_t j = 0 ; j < minC ; j++ )
1007  at( i,j ) = static_cast< T >( source_( i, j ));
1008 }
1009 
1010 
1011 template< size_t M, size_t N, typename T >
1012 void matrix< M, N, T >::operator=( const T old_fashioned_matrix[ M ][ N ] )
1013 {
1014  for( size_t row = 0; row < M; row++)
1015  {
1016  for( size_t col = 0; col < N; col++)
1017  {
1018  at( row, col ) = old_fashioned_matrix[ row ][ col ];
1019  }
1020  }
1021 
1022 }
1023 
1024 
1025 // WARNING: data_array[] must be at least of size M * N - otherwise CRASH!
1026 // WARNING: assumes row_by_row layout - if this is not the case,
1027 // use matrix::set( data_array, false )
1028 template< size_t M, size_t N, typename T >
1029 void
1030 matrix< M, N, T >::operator=( const T* data_array )
1031 {
1032  set( data_array, data_array + M * N, true );
1033 }
1034 
1035 
1036 
1037 template< size_t M, size_t N, typename T >
1038 void
1039 matrix< M, N, T >::operator=( const std::vector< T >& data )
1040 {
1041  if ( data.size() < M * N )
1042  throw std::runtime_error( "index out of bounds." );
1043 
1044  set( data.begin(), data.end(), true );
1045 }
1046 
1047 
1048 template< size_t M, size_t N, typename T >
1049 void
1050 matrix< M, N, T >::multiply_piecewise( const matrix& other )
1051 {
1052  for( size_t row_index = 0; row_index < M; row_index++)
1053  {
1054  for( size_t col_index = 0; col_index < N; col_index++ )
1055  {
1056  T& value = at( row_index, col_index );
1057  value *= other.at( row_index, col_index );
1058  }
1059  }
1060 }
1061 
1062 
1063 template< size_t M, size_t N, typename T >
1064 template< size_t P >
1065 void
1066 matrix< M, N, T >::multiply(
1067  const matrix< M, P, T >& left,
1068  const matrix< P, N, T >& right
1069  )
1070 {
1071  for( size_t row_index = 0; row_index < M; row_index++)
1072  {
1073  for( size_t col_index = 0; col_index < N; col_index++)
1074  {
1075  T& component = at( row_index, col_index );
1076  component = static_cast< T >( 0.0 );
1077  for( size_t p = 0; p < P; p++)
1078  {
1079  component += left.at( row_index, p ) * right.at( p, col_index );
1080  }
1081  }
1082  }
1083 }
1084 
1085 
1086 
1087 template< size_t M, size_t N, typename T >
1088 template< size_t P >
1089 matrix< M, P, T >
1090 matrix< M, N, T >::operator*( const matrix< N, P, T >& other ) const
1091 {
1092  matrix< M, P, T > result;
1093  result.multiply( *this, other );
1094  return result;
1095 }
1096 
1097 
1098 
1099 template< size_t M, size_t N, typename T >
1100 template< size_t O, size_t P, typename TT >
1101 typename enable_if< M == N && O == P && M == O, TT >::type*
1102 matrix< M, N, T >::operator*=( const matrix< O, P, TT >& right )
1103 {
1104  matrix< M, N, T > copy( *this );
1105  multiply( copy, right );
1106  return 0;
1107 }
1108 
1109 template< size_t M, size_t N, typename T >
1110 matrix< M, N, T >
1111 matrix< M, N, T >::operator/( T scalar )
1112 {
1113  matrix< M, N, T > result;
1114 
1115  for( size_t row_index = 0; row_index < M; ++row_index )
1116  {
1117  for( size_t col_index = 0; col_index < N; ++col_index )
1118  {
1119  result.at( row_index, col_index ) = at( row_index, col_index ) / scalar;
1120  }
1121  }
1122  return result;
1123 }
1124 
1125 
1126 
1127 template< size_t M, size_t N, typename T >
1128 void
1129 matrix< M, N, T >::operator/=( T scalar )
1130 {
1131  for( size_t row_index = 0; row_index < M; ++row_index )
1132  {
1133  for( size_t col_index = 0; col_index < N; ++col_index )
1134  {
1135  at( row_index, col_index ) /= scalar;
1136  }
1137  }
1138 }
1139 
1140 
1141 
1142 
1143 
1144 template< size_t M, size_t N, typename T >
1145 matrix< M, N, T >
1146 matrix< M, N, T >::operator*( T scalar )
1147 {
1148  matrix< M, N, T > result;
1149 
1150  for( size_t row_index = 0; row_index < M; ++row_index )
1151  {
1152  for( size_t col_index = 0; col_index < N; ++col_index )
1153  {
1154  result.at( row_index, col_index ) = at( row_index, col_index ) * scalar;
1155  }
1156  }
1157  return result;
1158 }
1159 
1160 
1161 
1162 template< size_t M, size_t N, typename T >
1163 void
1164 matrix< M, N, T >::operator*=( T scalar )
1165 {
1166  for( size_t row_index = 0; row_index < M; ++row_index )
1167  {
1168  for( size_t col_index = 0; col_index < N; ++col_index )
1169  {
1170  at( row_index, col_index ) *= scalar;
1171  }
1172  }
1173 }
1174 
1175 
1176 
1177 template< size_t M, size_t N, typename T >
1178 vector< M, T > matrix< M, N, T >::operator*( const vector< N, T >& vec ) const
1179 {
1180  vector< M, T > result;
1181 
1182  // this < M, 1 > = < M, P > * < P, 1 >
1183  for( size_t i = 0; i < M; ++i )
1184  {
1185  T tmp( 0 );
1186  for( size_t j = 0; j < N; ++j )
1187  tmp += at( i, j ) * vec.at( j );
1188  result.at( i ) = tmp;
1189  }
1190  return result;
1191 }
1192 
1193 
1194 
1195 // transform vector by matrix ( vec = matrix * vec )
1196 // assume homogenous coords( for vec3 = mat4 * vec3 ), e.g. vec[3] = 1.0
1197 template< size_t M, size_t N, typename T >
1198 template< size_t O > vector< O, T >
1199 matrix< M, N, T >::operator*( const vector< O, T >& vector_ ) const
1200 {
1201  vector< O, T > result;
1202  T tmp;
1203  for( size_t row_index = 0; row_index < M; ++row_index )
1204  {
1205  tmp = 0.0;
1206  for( size_t col_index = 0; col_index < N-1; ++col_index )
1207  {
1208  tmp += vector_( col_index ) * at( row_index, col_index );
1209  }
1210  if ( row_index < N - 1 )
1211  result( row_index ) = tmp + at( row_index, N-1 ); // * 1.0 -> homogeneous vec4
1212  else
1213  {
1214  tmp += at( row_index, N - 1 );
1215  for( size_t col_index = 0; col_index < N - 1; ++col_index )
1216  {
1217  result( col_index ) /= tmp;
1218  }
1219  }
1220  }
1221  return result;
1222 }
1223 
1224 
1225 
1226 template< size_t M, size_t N, typename T >
1227 inline matrix< M, N, T >
1228 matrix< M, N, T >::operator-() const
1229 {
1230  return negate();
1231 }
1232 
1233 
1234 
1235 template< size_t M, size_t N, typename T >
1236 matrix< M, N, T >
1237 matrix< M, N, T >::negate() const
1238 {
1239  matrix< M, N, T > result;
1240  result *= -1.0;
1241  return result;
1242 }
1243 
1244 
1245 
1246 template< size_t M, size_t N, typename T >
1247 void
1248 matrix< M, N, T >::tensor( const vector< M, T >& u, const vector< N, T >& v )
1249 {
1250  for ( size_t col_index = 0; col_index < N; ++col_index )
1251  for ( size_t row_index = 0; row_index < M; ++row_index )
1252  at( row_index, col_index ) = u.array[ col_index ] *
1253  v.array[ row_index ];
1254 }
1255 
1256 
1257 
1258 template< size_t M, size_t N, typename T >
1259 template< size_t uM, size_t vM >
1260 typename enable_if< uM == 3 && vM == 3 && M == N && M == 4 >::type*
1261 matrix< M, N, T >::tensor( const vector< uM, T >& u, const vector< vM, T >& v )
1262 {
1263  for ( size_t col_index = 0; col_index < 3; ++col_index )
1264  {
1265  for ( size_t row_index = 0; row_index < 3; ++row_index )
1266  at( row_index, col_index ) = u.array[ col_index ] *
1267  v.array[ row_index ];
1268 
1269  at( 3, col_index ) = u.array[ col_index ];
1270  }
1271 
1272  for ( size_t row_index = 0; row_index < 3; ++row_index )
1273  at( row_index, 3 ) = v.array[ row_index ];
1274 
1275  at( 3, 3 ) = 1.0;
1276 
1277  return 0;
1278 }
1279 
1280 
1281 
1282 template< size_t M, size_t N, typename T >
1283 void
1284 matrix< M, N, T >::
1285 transpose_to( matrix< N, M, T >& tM ) const
1286 {
1287  for( size_t row = 0; row < M; ++row )
1288  {
1289  for( size_t col = 0; col < N; ++col )
1290  {
1291  tM.at( col, row ) = at( row, col );
1292  }
1293  }
1294 }
1295 
1296 template< size_t M, size_t N, typename T >
1297 void
1298 matrix< M, N, T >::
1299 symmetric_covariance( matrix< M, M, T >& cov_m_ ) const
1300 {
1301  T tmp = 0;
1302  for( size_t row = 0; row < M; ++row )
1303  {
1304  for( size_t col = row; col < M; ++col )
1305  {
1306  for ( size_t k = 0; k < N; ++k )
1307  {
1308  tmp += (at( row, k ) * at( col, k ));
1309  }
1310 
1311  cov_m_.at( row, col ) = tmp;
1312  cov_m_.at( col, row ) = tmp;
1313  tmp = 0;
1314  }
1315  }
1316 }
1317 
1318 
1319 
1320 template< size_t M, size_t N, typename T >
1321 vector< M, T > matrix< M, N, T >::get_column( size_t index ) const
1322 {
1323  vector< M, T > column;
1324  get_column( index, column );
1325  return column;
1326 }
1327 
1328 
1329 
1330 template< size_t M, size_t N, typename T >
1331 void matrix< M, N, T >::get_column( size_t index, vector< M, T >& column ) const
1332 {
1333  if ( index >= N )
1334  throw std::runtime_error( "get_column() - index out of bounds." );
1335  memcpy( &column.array[0], &array[ M * index ], M * sizeof( T ) );
1336 }
1337 
1338 
1339 
1340 template< size_t M, size_t N, typename T >
1341 void matrix< M, N, T >::set_column( size_t index, const vector< M, T >& column )
1342 {
1343  if ( index >= N )
1344  throw std::runtime_error( "set_column() - index out of bounds." );
1345  memcpy( array + M * index, column.array, M * sizeof( T ) );
1346 }
1347 
1348 template< size_t M, size_t N, typename T >
1349 void matrix< M, N, T >::get_column( size_t index, matrix< M, 1, T >& column )
1350  const
1351 {
1352  if ( index >= N )
1353  throw std::runtime_error( "get_column() - index out of bounds." );
1354  memcpy( column.array, array + M * index, M * sizeof( T ) );
1355 }
1356 
1357 template< size_t M, size_t N, typename T >
1358 void matrix< M, N, T >::set_column( size_t index,
1359  const matrix< M, 1, T >& column )
1360 {
1361  if ( index >= N )
1362  throw std::runtime_error( "set_column() - index out of bounds." );
1363  memcpy( &array[ M * index ], column.array, M * sizeof( T ) );
1364 }
1365 
1366 template< size_t M, size_t N, typename T >
1367 vector< N, T > matrix< M, N, T >::get_row( size_t index ) const
1368 {
1369  vector< N, T > row;
1370  get_row( index, row );
1371  return row;
1372 }
1373 
1374 template< size_t M, size_t N, typename T >
1375 void matrix< M, N, T >::get_row( size_t row_index, vector< N, T >& row ) const
1376 {
1377  if ( row_index >= M )
1378  throw std::runtime_error( "get_row() - index out of bounds." );
1379 
1380  for( size_t col_index = 0; col_index < N; ++col_index )
1381  row.at( col_index ) = at( row_index, col_index );
1382 }
1383 
1384 template< size_t M, size_t N, typename T >
1385 void matrix< M, N, T >::set_row( size_t row_index, const vector< N, T >& row )
1386 {
1387  if ( row_index >= M )
1388  throw std::runtime_error( "set_row() - index out of bounds." );
1389 
1390  for( size_t col_index = 0; col_index < N; ++col_index )
1391  at( row_index, col_index ) = row.at( col_index );
1392 }
1393 
1394 template< size_t M, size_t N, typename T >
1395 void matrix< M, N, T >::get_row( size_t row_index, matrix< 1, N, T >& row )
1396  const
1397 {
1398  if ( row_index >= M )
1399  throw std::runtime_error( "get_row() - index out of bounds." );
1400 
1401  for( size_t col_index = 0; col_index < N; ++col_index )
1402  row.at( 0, col_index ) = at( row_index, col_index );
1403 }
1404 
1405 template< size_t M, size_t N, typename T >
1406 void matrix< M, N, T >::set_row( size_t row_index,
1407  const matrix< 1, N, T >& row )
1408 {
1409  if ( row_index >= M )
1410  throw std::runtime_error( "set_row() - index out of bounds." );
1411 
1412  for( size_t col_index = 0; col_index < N; ++col_index )
1413  at( row_index, col_index ) = row.at( 0, col_index );
1414 }
1415 
1416 template< size_t M, size_t N, typename T >
1417 size_t matrix< M, N, T >::get_number_of_rows() const
1418 {
1419  return M;
1420 }
1421 
1422 template< size_t M, size_t N, typename T >
1423 size_t matrix< M, N, T >::get_number_of_columns() const
1424 {
1425  return N;
1426 }
1427 
1428 template< size_t M, size_t N, typename T >
1429 void
1430 matrix< M, N, T >::
1431 fill( T fillValue )
1432 {
1433  for( size_t row_index = 0; row_index < M; ++row_index )
1434  {
1435  for( size_t col_index = 0; col_index < N; ++col_index )
1436  {
1437  at( row_index, col_index ) = fillValue;
1438  }
1439  }
1440 }
1441 
1442 
1443 template< size_t M, size_t N, typename T >
1444 inline T&
1445 matrix< M, N, T >::x()
1446 {
1447  return array[ 12 ];
1448 }
1449 
1450 
1451 
1452 template< size_t M, size_t N, typename T >
1453 inline T&
1454 matrix< M, N, T >::y()
1455 {
1456  return array[ 13 ];
1457 }
1458 
1459 
1460 
1461 template< size_t M, size_t N, typename T >
1462 inline T&
1463 matrix< M, N, T >::z()
1464 {
1465  return array[ 14 ];
1466 }
1467 
1468 template< size_t M, size_t N, typename T >
1469 template< typename input_iterator_t >
1470 void matrix< M, N, T >::set( input_iterator_t begin_,
1471  input_iterator_t end_, bool row_major_layout )
1472 {
1473  input_iterator_t it( begin_ );
1474  if( row_major_layout )
1475  {
1476  for( size_t row = 0; row < M; ++row )
1477  {
1478  for( size_t col = 0; col < N; ++col, ++it )
1479  {
1480  if ( it == end_ )
1481  return;
1482  at( row, col ) = static_cast< T >( *it );
1483  }
1484  }
1485  }
1486  else
1487  {
1488  std::copy( it, it + ( M * N ), begin() );
1489  }
1490 }
1491 
1492 
1493 
1494 template< size_t M, size_t N, typename T >
1495 void matrix< M, N, T >::zero()
1496 {
1497  fill( static_cast< T >( 0.0 ) );
1498 }
1499 
1500 template< size_t M, size_t N, typename T >
1501 void matrix< M, N, T >::operator=( T value_ )
1502 {
1503  std::fill( begin(), end(), value_ );
1504 }
1505 
1506 
1507 
1508 template< size_t M, size_t N, typename T >
1509 inline matrix< M, N, T >
1510 matrix< M, N, T >::operator+( const matrix< M, N, T >& other ) const
1511 {
1512  matrix< M, N, T > result( *this );
1513  result += other;
1514  return result;
1515 }
1516 
1517 
1518 
1519 template< size_t M, size_t N, typename T >
1520 void matrix< M, N, T >::operator+=( const matrix< M, N, T >& other )
1521 {
1522  iterator it = begin(), it_end = end();
1523  const_iterator other_it = other.begin();
1524  for( ; it != it_end; ++it, ++other_it )
1525  {
1526  *it += *other_it;
1527  }
1528 }
1529 
1530 
1531 template< size_t M, size_t N, typename T >
1532 void matrix< M, N, T >::operator+=( T scalar )
1533 {
1534  iterator it = begin(), it_end = end();
1535  for( ; it != it_end; ++it )
1536  {
1537  *it += scalar;
1538  }
1539 }
1540 
1541 template< size_t M, size_t N, typename T >
1542 inline matrix< M, N, T >
1543 matrix< M, N, T >::
1544 operator-( const matrix< M, N, T >& other ) const
1545 {
1546  matrix< M, N, T > result( *this );
1547  result -= other;
1548  return result;
1549 }
1550 
1551 
1552 
1553 template< size_t M, size_t N, typename T >
1554 void
1555 matrix< M, N, T >::
1556 operator-=( const matrix< M, N, T >& other )
1557 {
1558  iterator it = begin(), it_end = end();
1559  const_iterator other_it = other.begin();
1560  for( ; it != it_end; ++it, ++other_it )
1561  {
1562  *it -= *other_it;
1563  }
1564 }
1565 
1566 template< size_t M, size_t N, typename T >
1567 void
1568 matrix< M, N, T >::operator-=( T scalar )
1569 {
1570  iterator it = begin(), it_end = end();
1571  for( ; it != it_end; ++it )
1572  {
1573  *it -= scalar;
1574  }
1575 }
1576 
1577 
1578 template< size_t M, size_t N, typename T >
1579 template< size_t O, size_t P, size_t Q, size_t R >
1580 typename enable_if< M == O + Q && N == P + R >::type*
1581 matrix< M, N, T >::direct_sum( const matrix< O, P, T >& upper_left,
1582  const matrix< Q, R, T >& lower_right )
1583 {
1584  (*this) = static_cast< T >( 0.0 );
1585 
1586  for( size_t row = 0; row < O; ++row )
1587  {
1588  for( size_t col = 0; col < P; ++col )
1589  {
1590  at( row, col ) = upper_left( row, col );
1591  }
1592  }
1593 
1594  for( size_t row = 0; row < Q; ++row )
1595  {
1596  for( size_t col = 0; col < R; ++col )
1597  {
1598  at( O + row, P + col ) = lower_right( row, col );
1599  }
1600  }
1601 
1602  return 0;
1603 }
1604 
1605 
1606 
1607 template< size_t M, size_t N, typename T >
1608 template< size_t O, size_t P >
1609 matrix< O, P, T >
1610 matrix< M, N, T >::get_sub_matrix( size_t row_offset, size_t col_offset,
1611 typename enable_if< O <= M && P <= N >::type* ) const
1612 {
1613  matrix< O, P, T > result;
1614  get_sub_matrix( result, row_offset, col_offset );
1615  return result;
1616 }
1617 
1618 
1619 
1620 template< size_t M, size_t N, typename T >
1621 template< size_t O, size_t P >
1622 typename enable_if< O <= M && P <= N >::type*
1623 matrix< M, N, T >::
1624 get_sub_matrix( matrix< O, P, T >& result,
1625  size_t row_offset, size_t col_offset ) const
1626 {
1627  #ifdef VMMLIB_SAFE_ACCESSORS
1628  if ( O + row_offset > M || P + col_offset > N )
1629  throw std::runtime_error( "index out of bounds." );
1630  #endif
1631 
1632  for( size_t row = 0; row < O; ++row )
1633  {
1634  for( size_t col = 0; col < P; ++col )
1635  {
1636  result.at( row, col )
1637  = at( row_offset + row, col_offset + col );
1638  }
1639  }
1640  return 0;
1641 }
1642 
1643 
1644 
1645 template< size_t M, size_t N, typename T >
1646 template< size_t O, size_t P >
1647 typename enable_if< O <= M && P <= N >::type*
1648 matrix< M, N, T >::
1649 set_sub_matrix( const matrix< O, P, T >& sub_matrix,
1650  size_t row_offset, size_t col_offset )
1651 {
1652  for( size_t row = 0; row < O; ++row )
1653  {
1654  for( size_t col = 0; col < P; ++col )
1655  {
1656  at( row_offset + row, col_offset + col )
1657  = sub_matrix.at( row, col );
1658  }
1659  }
1660  return 0; // for sfinae
1661 }
1662 
1663 
1664 
1665 template< size_t M, size_t N, typename T >
1666 inline T
1667 matrix< M, N, T >::det() const
1668 {
1669  return compute_determinant( *this );
1670 }
1671 
1672 
1673 
1674 template< size_t M, size_t N, typename T >
1675 template< size_t O, size_t P, typename TT >
1676 inline bool matrix< M, N, T >::inverse( matrix< O, P, TT >& inverse_,
1677  T tolerance, typename
1678  enable_if< M == N && O == P && O == M && M >= 2 && M <= 4, TT >::type* )
1679  const
1680 {
1681  return compute_inverse( *this, inverse_, tolerance );
1682 }
1683 
1684 
1685 
1686 template< size_t M, size_t N, typename T >
1687 template< size_t O, size_t P >
1688 typename enable_if< O == P && M == N && O == M && M >= 2 >::type*
1689 matrix< M, N, T >::
1690 get_adjugate( matrix< O, P, T >& adjugate ) const
1691 {
1692  get_cofactors( adjugate );
1693  adjugate = transpose( adjugate );
1694  return 0;
1695 }
1696 
1697 
1698 
1699 template< size_t M, size_t N, typename T >
1700 template< size_t O, size_t P >
1701 typename enable_if< O == P && M == N && O == M && M >= 2 >::type*
1702 matrix< M, N, T >::
1703 get_cofactors( matrix< O, P, T >& cofactors ) const
1704 {
1705  matrix< M-1, N-1, T > minor_;
1706 
1707  const size_t _negate = 1u;
1708  for( size_t row_index = 0; row_index < M; ++row_index )
1709  {
1710  for( size_t col_index = 0; col_index < N; ++col_index )
1711  {
1712  if ( ( row_index + col_index ) & _negate )
1713  cofactors( row_index, col_index ) = -get_minor( minor_, row_index, col_index );
1714  else
1715  cofactors( row_index, col_index ) = get_minor( minor_, row_index, col_index );
1716  }
1717  }
1718  return 0;
1719 }
1720 
1721 
1722 
1723 template< size_t M, size_t N, typename T >
1724 template< size_t O, size_t P >
1725 T
1726 matrix< M, N, T >::
1727 get_minor( matrix< O, P, T >& minor_, size_t row_to_cut, size_t col_to_cut,
1728 typename enable_if< O == M-1 && P == N-1 && M == N && M >= 2 >::type* ) const
1729 {
1730  size_t row_offset = 0;
1731  size_t col_offset = 0;
1732  for( size_t row_index = 0; row_index < M; ++row_index )
1733  {
1734  if ( row_index == row_to_cut )
1735  row_offset = -1;
1736  else
1737  {
1738  for( size_t col_index = 0; col_index < M; ++col_index )
1739  {
1740  if ( col_index == col_to_cut )
1741  col_offset = -1;
1742  else
1743  minor_.at( row_index + row_offset, col_index + col_offset )
1744  = at( row_index, col_index );
1745  }
1746  col_offset = 0;
1747  }
1748  }
1749  return compute_determinant( minor_ );
1750 }
1751 
1752 template< size_t M, size_t N, typename T >
1753 template< typename TT >
1754 matrix< M, N, T >& matrix< M, N, T >::rotate(
1755  const TT angle_, const vector< M-1, T >& axis,
1756  typename enable_if< M == N && M == 4, TT >::type* )
1757 {
1758  const T angle = static_cast< T >( angle_ );
1759 
1760  const T sine = std::sin( angle );
1761  const T cosine = std::cos( angle );
1762 
1763  // this is necessary since Visual Studio cannot resolve the
1764  // std::pow()-call correctly if we just use 2.0 directly.
1765  // this way, the '2.0' is converted to the same format
1766  // as the axis components
1767 
1768  const T _zero = 0.0;
1769  const T one = 1.0;
1770  const T two = 2.0;
1771 
1772  array[0] = cosine + ( one - cosine ) * std::pow( axis.array[0], two );
1773  array[1] = ( one - cosine ) * axis.array[0] * axis.array[1] + sine * axis.array[2];
1774  array[2] = ( one - cosine ) * axis.array[0] * axis.array[2] - sine * axis.array[1];
1775  array[3] = _zero;
1776 
1777  array[4] = ( one - cosine ) * axis.array[0] * axis.array[1] - sine * axis.array[2];
1778  array[5] = cosine + ( one - cosine ) * std::pow( axis.array[1], two );
1779  array[6] = ( one - cosine ) * axis.array[1] * axis.array[2] + sine * axis.array[0];
1780  array[7] = _zero;
1781 
1782  array[8] = ( one - cosine ) * axis.array[0] * axis.array[2] + sine * axis.array[1];
1783  array[9] = ( one - cosine ) * axis.array[1] * axis.array[2] - sine * axis.array[0];
1784  array[10] = cosine + ( one - cosine ) * std::pow( axis.array[2], two );
1785  array[11] = _zero;
1786 
1787  array[12] = _zero;
1788  array[13] = _zero;
1789  array[14] = _zero;
1790  array[15] = one;
1791 
1792  return *this;
1793 }
1794 
1795 template< size_t M, size_t N, typename T >
1796 template< typename TT >
1797 matrix< M, N, T >& matrix< M, N, T >::rotate_x( const TT angle_,
1798  typename enable_if< M == N && M == 4, TT >::type* )
1799 {
1800  const T angle = static_cast< T >( angle_ );
1801 
1802  const T sine = std::sin( angle );
1803  const T cosine = std::cos( angle );
1804 
1805  T tmp;
1806 
1807  tmp = array[ 4 ] * cosine + array[ 8 ] * sine;
1808  array[ 8 ] = - array[ 4 ] * sine + array[ 8 ] * cosine;
1809  array[ 4 ] = tmp;
1810 
1811  tmp = array[ 5 ] * cosine + array[ 9 ] * sine;
1812  array[ 9 ] = - array[ 5 ] * sine + array[ 9 ] * cosine;
1813  array[ 5 ] = tmp;
1814 
1815  tmp = array[ 6 ] * cosine + array[ 10 ] * sine;
1816  array[ 10 ] = - array[ 6 ] * sine + array[ 10 ] * cosine;
1817  array[ 6 ] = tmp;
1818 
1819  tmp = array[ 7 ] * cosine + array[ 11 ] * sine;
1820  array[ 11 ] = - array[ 7 ] * sine + array[ 11 ] * cosine;
1821  array[ 7 ] = tmp;
1822 
1823  return *this;
1824 }
1825 
1826 template< size_t M, size_t N, typename T >
1827 template< typename TT >
1828 matrix< M, N, T >& matrix< M, N, T >::rotate_y( const TT angle_,
1829  typename enable_if< M == N && M == 4, TT >::type* )
1830 {
1831  const T angle = static_cast< T >( angle_ );
1832 
1833  const T sine = std::sin( angle );
1834  const T cosine = std::cos( angle );
1835 
1836  T tmp;
1837 
1838  tmp = array[ 0 ] * cosine - array[ 8 ] * sine;
1839  array[ 8 ] = array[ 0 ] * sine + array[ 8 ] * cosine;
1840  array[ 0 ] = tmp;
1841 
1842  tmp = array[ 1 ] * cosine - array[ 9 ] * sine;
1843  array[ 9 ] = array[ 1 ] * sine + array[ 9 ] * cosine;
1844  array[ 1 ] = tmp;
1845 
1846  tmp = array[ 2 ] * cosine - array[ 10 ] * sine;
1847  array[ 10 ] = array[ 2 ] * sine + array[ 10 ] * cosine;
1848  array[ 2 ] = tmp;
1849 
1850  tmp = array[ 3 ] * cosine - array[ 11 ] * sine;
1851  array[ 11 ] = array[ 3 ] * sine + array[ 11 ] * cosine;
1852  array[ 3 ] = tmp;
1853 
1854  return *this;
1855 }
1856 
1857 template< size_t M, size_t N, typename T >
1858 template< typename TT >
1859 matrix< M, N, T >& matrix< M, N, T >::rotate_z( const TT angle_,
1860  typename enable_if< M == N && M == 4, TT >::type* )
1861 {
1862  const T angle = static_cast< T >( angle_ );
1863 
1864  const T sine = std::sin( angle );
1865  const T cosine = std::cos( angle );
1866 
1867  T tmp;
1868 
1869  tmp = array[ 0 ] * cosine + array[ 4 ] * sine;
1870  array[ 4 ] = - array[ 0 ] * sine + array[ 4 ] * cosine;
1871  array[ 0 ] = tmp;
1872 
1873  tmp = array[ 1 ] * cosine + array[ 5 ] * sine;
1874  array[ 5 ] = - array[ 1 ] * sine + array[ 5 ] * cosine;
1875  array[ 1 ] = tmp;
1876 
1877  tmp = array[ 2 ] * cosine + array[ 6 ] * sine;
1878  array[ 6 ] = - array[ 2 ] * sine + array[ 6 ] * cosine;
1879  array[ 2 ] = tmp;
1880 
1881  tmp = array[ 3 ] * cosine + array[ 7 ] * sine;
1882  array[ 7 ] = - array[ 3 ] * sine + array[ 7 ] * cosine;
1883  array[ 3 ] = tmp;
1884 
1885  return *this;
1886 }
1887 
1888 template< size_t M, size_t N, typename T >
1889 template< typename TT >
1890 matrix< M, N, T >& matrix< M, N, T >::pre_rotate_x( const TT angle_,
1891  typename enable_if< M == N && M == 4, TT >::type* )
1892 {
1893  const T angle = static_cast< T >( angle_ );
1894 
1895  const T sine = std::sin( angle );
1896  const T cosine = std::cos( angle );
1897 
1898  T tmp;
1899 
1900  tmp = array[ 1 ];
1901  array[ 1 ] = array[ 1 ] * cosine + array[ 2 ] * sine;
1902  array[ 2 ] = tmp * -sine + array[ 2 ] * cosine;
1903 
1904  tmp = array[ 5 ];
1905  array[ 5 ] = array[ 5 ] * cosine + array[ 6 ] * sine;
1906  array[ 6 ] = tmp * -sine + array[ 6 ] * cosine;
1907 
1908  tmp = array[ 9 ];
1909  array[ 9 ] = array[ 9 ] * cosine + array[ 10 ] * sine;
1910  array[ 10 ] = tmp * -sine + array[ 10 ] * cosine;
1911 
1912  tmp = array[ 13 ];
1913  array[ 13 ] = array[ 13 ] * cosine + array[ 14 ] * sine;
1914  array[ 14 ] = tmp * -sine + array[ 14 ] * cosine;
1915 
1916  return *this;
1917 }
1918 
1919 template< size_t M, size_t N, typename T >
1920 template< typename TT >
1921 matrix< M, N, T >& matrix< M, N, T >::pre_rotate_y( const TT angle_,
1922  typename enable_if< M == N && M == 4, TT >::type* )
1923 {
1924  const T angle = static_cast< T >( angle_ );
1925 
1926  const T sine = std::sin( angle );
1927  const T cosine = std::cos( angle );
1928 
1929  T tmp;
1930 
1931  tmp = array[ 0 ];
1932  array[ 0 ] = array[ 0 ] * cosine - array[ 2 ] * sine;
1933  array[ 2 ] = tmp * sine + array[ 2 ] * cosine;
1934 
1935  tmp = array[ 4 ];
1936  array[ 4 ] = array[ 4 ] * cosine - array[ 6 ] * sine;
1937  array[ 6 ] = tmp * sine + array[ 6 ] * cosine;
1938 
1939  tmp = array[ 8 ];
1940  array[ 8 ] = array[ 8 ] * cosine - array[ 10 ] * sine;
1941  array[ 10 ] = tmp * sine + array[ 10 ] * cosine;
1942 
1943  tmp = array[ 12 ];
1944  array[ 12 ] = array[ 12 ] * cosine - array[ 14 ] * sine;
1945  array[ 14 ] = tmp * sine + array[ 14 ] * cosine;
1946 
1947  return *this;
1948 }
1949 
1950 template< size_t M, size_t N, typename T >
1951 template< typename TT >
1952 matrix< M, N, T >& matrix< M, N, T >::pre_rotate_z( const TT angle_,
1953  typename enable_if< M == N && M == 4, TT >::type* )
1954 {
1955  const T angle = static_cast< T >( angle_ );
1956 
1957  const T sine = std::sin( angle );
1958  const T cosine = std::cos( angle );
1959 
1960  T tmp;
1961 
1962  tmp = array[ 0 ];
1963  array[ 0 ] = array[ 0 ] * cosine + array[ 1 ] * sine;
1964  array[ 1 ] = tmp * -sine + array[ 1 ] * cosine;
1965 
1966  tmp = array[ 4 ];
1967  array[ 4 ] = array[ 4 ] * cosine + array[ 5 ] * sine;
1968  array[ 5 ] = tmp * -sine + array[ 5 ] * cosine;
1969 
1970  tmp = array[ 8 ];
1971  array[ 8 ] = array[ 8 ] * cosine + array[ 9 ] * sine;
1972  array[ 9 ] = tmp * -sine + array[ 9 ] * cosine;
1973 
1974  tmp = array[ 12 ];
1975  array[ 12 ] = array[ 12 ] * cosine + array[ 13 ] * sine;
1976  array[ 13 ] = tmp * -sine + array[ 13 ] * cosine;
1977 
1978  return *this;
1979 }
1980 
1981 template< size_t M, size_t N, typename T >
1982 template< typename TT >
1983 matrix< M, N, T >& matrix< M, N, T >::scale( const TT _scale[3],
1984  typename enable_if< M == N && M == 4, TT >::type* )
1985 {
1986  const T scale0 = static_cast< T >( _scale[ 0 ] );
1987  const T scale1 = static_cast< T >( _scale[ 1 ] );
1988  const T scale2 = static_cast< T >( _scale[ 2 ] );
1989 
1990  array[0] *= scale0;
1991  array[1] *= scale0;
1992  array[2] *= scale0;
1993  array[3] *= scale0;
1994  array[4] *= scale1;
1995  array[5] *= scale1;
1996  array[6] *= scale1;
1997  array[7] *= scale1;
1998  array[8] *= scale2;
1999  array[9] *= scale2;
2000  array[10] *= scale2;
2001  array[11] *= scale2;
2002 
2003  return *this;
2004 }
2005 
2006 template< size_t M, size_t N, typename T >
2007 template< typename TT >
2008 matrix< M, N, T >& matrix< M, N, T >::scale( const TT x_, const T y_, const T z_,
2009  typename enable_if< M == N && M == 4, TT >::type* )
2010 {
2011  const T _x = static_cast< T >( x_ );
2012 
2013  array[0] *= _x;
2014  array[1] *= _x;
2015  array[2] *= _x;
2016  array[3] *= _x;
2017  array[4] *= y_;
2018  array[5] *= y_;
2019  array[6] *= y_;
2020  array[7] *= y_;
2021  array[8] *= z_;
2022  array[9] *= z_;
2023  array[10] *= z_;
2024  array[11] *= z_;
2025 
2026  return *this;
2027 }
2028 
2029 template< size_t M, size_t N, typename T >
2030 template< typename TT >
2031 inline matrix< M, N, T >& matrix< M, N, T >::scale(
2032  const vector< 3, TT >& scale_,
2033  typename enable_if< M == N && M == 4, TT >::type* )
2034 {
2035  return scale( scale_.array );
2036 }
2037 
2038 
2039 
2040 template< size_t M, size_t N, typename T >
2041 template< typename TT >
2042 matrix< M, N, T >& matrix< M, N, T >::scale_translation( const TT scale_[3],
2043  typename enable_if< M == N && M == 4, TT >::type* )
2044 {
2045  array[12] *= static_cast< T >( scale_[0] );
2046  array[13] *= static_cast< T >( scale_[1] );
2047  array[14] *= static_cast< T >( scale_[2] );
2048  return *this;
2049 }
2050 
2051 template< size_t M, size_t N, typename T >
2052 template< typename TT >
2053 inline matrix< M, N, T >& matrix< M, N, T >::scale_translation(
2054  const vector< 3, TT >& scale_,
2055  typename enable_if< M == N && M == 4, TT >::type* )
2056 {
2057  return scale_translation( scale_.array );
2058 }
2059 
2060 template< size_t M, size_t N, typename T >
2061 template< typename TT >
2062 inline matrix< M, N, T >& matrix< M, N, T >::set_translation(
2063  const TT x_, const TT y_, const TT z_,
2064  typename enable_if< M == N && M == 4, TT >::type* )
2065 {
2066  array[12] = static_cast< T >( x_ );
2067  array[13] = static_cast< T >( y_ );
2068  array[14] = static_cast< T >( z_ );
2069  return *this;
2070 }
2071 
2072 template< size_t M, size_t N, typename T >
2073 template< typename TT >
2074 inline matrix< M, N, T >& matrix< M, N, T >::set_translation( const TT trans[3],
2075  typename enable_if< M == N && M == 4, TT >::type* )
2076 {
2077  array[12] = static_cast< T >( trans[ 0 ] );
2078  array[13] = static_cast< T >( trans[ 1 ] );
2079  array[14] = static_cast< T >( trans[ 2 ] );
2080  return *this;
2081 }
2082 
2083 template< size_t M, size_t N, typename T >
2084 template< typename TT >
2085 inline matrix< M, N, T >& matrix< M, N, T >::set_translation(
2086  const vector< 3, TT >& translation_,
2087  typename enable_if< M == N && M == 4, TT >::type* )
2088 {
2089  return set_translation( translation_.array );
2090 }
2091 
2092 template< size_t M, size_t N, typename T > template< typename TT > inline
2093 void matrix< M, N, T >::get_translation( vector< N-1, TT >& translation_ )
2094  const
2095 {
2096  for( size_t i = 0; i < N-1; ++i )
2097  translation_.array[ i ] = array[ i + M * (N - 1) ];
2098 }
2099 
2100 
2101 template< size_t M, size_t N, typename T >
2102 inline vector< N-1, T > matrix< M, N, T >::get_translation() const
2103 {
2104  vector< N-1, T > result;
2105  get_translation( result );
2106  return result;
2107 }
2108 
2109 template< size_t M, size_t N, typename T >
2110 template< size_t O >
2111 void matrix< M, N, T >::getLookAt( vector< O, T >& eye, vector< O, T >& lookAt,
2112  vector< O, T >& up,
2113  typename enable_if< M == O+1 && N == O+1 && O == 3 >::type* ) const
2114 {
2115  matrix< 4, 4, T > inv;
2116  inverse( inv );
2117 
2118  const matrix< 3, 3, T > rotation( transpose( matrix< 3, 3, T >( *this )));
2119 
2120  eye = inv * vector< 3, T >::ZERO;
2121  up = rotation * vector< 3, T >::UP;
2122  lookAt = rotation * vector< 3, T >::FORWARD;
2123  lookAt.normalize();
2124  lookAt = eye + lookAt;
2125 }
2126 
2127 template< size_t M, size_t N, typename T >
2128 size_t
2129 matrix< M, N, T >::
2130 size() const
2131 {
2132  return M * N;
2133 }
2134 
2135 
2136 
2137 template< size_t M, size_t N, typename T >
2138 typename matrix< M, N, T >::iterator
2139 matrix< M, N, T >::
2140 begin()
2141 {
2142  return array;
2143 }
2144 
2145 
2146 
2147 template< size_t M, size_t N, typename T >
2148 typename matrix< M, N, T >::iterator
2149 matrix< M, N, T >::
2150 end()
2151 {
2152  return array + size();
2153 }
2154 
2155 
2156 
2157 template< size_t M, size_t N, typename T >
2158 typename matrix< M, N, T >::const_iterator
2159 matrix< M, N, T >::
2160 begin() const
2161 {
2162  return array;
2163 }
2164 
2165 
2166 
2167 template< size_t M, size_t N, typename T >
2168 typename matrix< M, N, T >::const_iterator
2169 matrix< M, N, T >::
2170 end() const
2171 {
2172  return array + size();
2173 }
2174 
2175 
2176 
2177 template< size_t M, size_t N, typename T >
2178 typename matrix< M, N, T >::reverse_iterator
2179 matrix< M, N, T >::
2180 rbegin()
2181 {
2182  return array + size() - 1;
2183 }
2184 
2185 
2186 
2187 template< size_t M, size_t N, typename T >
2188 typename matrix< M, N, T >::reverse_iterator
2189 matrix< M, N, T >::
2190 rend()
2191 {
2192  return array - 1;
2193 }
2194 
2195 
2196 
2197 template< size_t M, size_t N, typename T >
2198 typename matrix< M, N, T >::const_reverse_iterator
2199 matrix< M, N, T >::
2200 rbegin() const
2201 {
2202  return array + size() - 1;
2203 }
2204 
2205 
2206 
2207 template< size_t M, size_t N, typename T >
2208 typename matrix< M, N, T >::const_reverse_iterator
2209 matrix< M, N, T >::
2210 rend() const
2211 {
2212  return array - 1;
2213 }
2214 
2215 
2216 
2217 template< size_t M, size_t N, typename T > template< typename init_functor_t >
2218 const matrix< M, N, T > matrix< M, N, T >::get_initialized_matrix()
2219 {
2220  matrix< M, N, T > matrix_;
2221  init_functor_t()( matrix_ );
2222  return matrix_;
2223 }
2224 
2225 
2226 // it's ugly, but it allows having properly initialized static members
2227 // without having to worry about static initialization order.
2228 template< size_t M, size_t N, typename T >
2229 const matrix< M, N, T >
2230 matrix< M, N, T >::IDENTITY(
2231  matrix< M, N, T >::get_initialized_matrix<
2232  set_to_identity_functor< matrix< M, N, T > > >( ));
2233 
2234 
2235 template< size_t M, size_t N, typename T >
2236 const matrix< M, N, T >
2237 matrix< M, N, T >::ZERO(
2238  matrix< M, N, T >::
2239  get_initialized_matrix< set_to_zero_functor< matrix< M, N, T > > >()
2240  );
2241 
2242 
2243 
2244 template< size_t M, size_t N, typename T >
2245 double
2246 matrix< M, N, T >::frobenius_norm( ) const
2247 {
2248  double norm = 0.0;
2249 
2250  const_iterator it = begin(), it_end = end();
2251  for( ; it != it_end; ++it )
2252  {
2253  norm += *it * *it;
2254  }
2255 
2256  return std::sqrt(norm);
2257 }
2258 
2259 template< size_t M, size_t N, typename T >
2260 double
2261 matrix< M, N, T >::p_norm( double p ) const
2262 {
2263  double norm = 0.0;
2264 
2265  const_iterator it = begin(), it_end = end();
2266  for( ; it != it_end; ++it )
2267  {
2268  norm += std::pow(*it, p);
2269  }
2270 
2271  return std::pow(norm,1./p);
2272 }
2273 
2274 
2275 template< size_t M, size_t N, typename T >
2276 template< size_t O >
2277 void
2278 matrix< M, N, T >::khatri_rao_product( const matrix< O, N, T >& right_, matrix< M*O, N, T >& prod_ ) const
2279 {
2280  //build product for every column
2281  for (size_t col = 0; col < N; ++col )
2282  {
2283  for ( size_t m = 0; m < M; ++m )
2284  {
2285  for (size_t o = 0; o < O; ++o )
2286  {
2287  prod_.at(O*m + o, col) = at( m, col ) * right_.at( o, col );
2288  }
2289  }
2290  }
2291 }
2292 
2293 template< size_t M, size_t N, typename T >
2294 template< size_t O, size_t P >
2295 void
2296 matrix< M, N, T >::kronecker_product( const matrix< O, P, T >& right_, matrix< M*O, N*P, T >& result_ ) const
2297 {
2298  //build product for every column
2299  for (size_t m = 0; m < M; ++m )
2300  {
2301  for ( size_t n = 0; n < N; ++n )
2302  {
2303  for (size_t o = 0; o < O; ++o )
2304  {
2305  for (size_t p = 0; p < P; ++p )
2306  {
2307  result_.at(O*m + o, P*n + p) = at( m, n ) * right_.at( o, p );
2308  }
2309  }
2310  }
2311  }
2312 }
2313 
2314 
2315 template< size_t M, size_t N, typename T >
2316 template< typename TT >
2317 void
2318 matrix< M, N, T >::cast_from( const matrix< M, N, TT >& other )
2319 {
2320  typedef vmml::matrix< M, N, TT > matrix_tt_type ;
2321  typedef typename matrix_tt_type::const_iterator tt_const_iterator;
2322 
2323  iterator it = begin(), it_end = end();
2324  tt_const_iterator other_it = other.begin();
2325  for( ; it != it_end; ++it, ++other_it )
2326  {
2327  *it = static_cast< T >( *other_it );
2328  }
2329 }
2330 
2331 
2332 template< size_t M, size_t N, typename T >
2333 T
2334 matrix< M, N, T >::get_min() const
2335 {
2336  T min_value = static_cast<T>(std::numeric_limits<T>::max());
2337 
2338  const_iterator it = begin(),
2339  it_end = end();
2340  for( ; it != it_end; ++it)
2341  {
2342  if ( *it < min_value ) {
2343  min_value = *it;
2344  }
2345  }
2346  return min_value;
2347 }
2348 
2349 template< size_t M, size_t N, typename T >
2350 T
2351 matrix< M, N, T >::get_max() const
2352 {
2353  T max_value = static_cast<T>(0);
2354 
2355  const_iterator it = begin(),
2356  it_end = end();
2357  for( ; it != it_end; ++it)
2358  {
2359  if ( *it > max_value ) {
2360  max_value = *it;
2361  }
2362  }
2363  return max_value;
2364 }
2365 
2366 
2367 template< size_t M, size_t N, typename T >
2368 T
2369 matrix< M, N, T >::get_abs_min() const
2370 {
2371  T min_value = static_cast<T>(std::numeric_limits<T>::max());
2372 
2373  const_iterator it = begin(),
2374  it_end = end();
2375  for( ; it != it_end; ++it)
2376  {
2377  if ( fabs(*it) < fabs(min_value) ) {
2378  min_value = fabs(*it);
2379  }
2380  }
2381  return min_value;
2382 }
2383 
2384 template< size_t M, size_t N, typename T >
2385 T
2386 matrix< M, N, T >::get_abs_max() const
2387 {
2388  T max_value = static_cast<T>(0);
2389 
2390  const_iterator it = begin(),
2391  it_end = end();
2392  for( ; it != it_end; ++it)
2393  {
2394  if ( fabs(*it) > fabs(max_value) ) {
2395  max_value = fabs(*it);
2396  }
2397  }
2398  return max_value;
2399 }
2400 
2401 
2402 
2403 template< size_t M, size_t N, typename T >
2404 size_t
2405 matrix< M, N, T >::nnz() const
2406 {
2407  size_t counter = 0;
2408 
2409  const_iterator it = begin(),
2410  it_end = end();
2411  for( ; it != it_end; ++it)
2412  {
2413  if ( *it != 0 ) {
2414  ++counter;
2415  }
2416  }
2417 
2418  return counter;
2419 }
2420 
2421 template< size_t M, size_t N, typename T >
2422 size_t
2423 matrix< M, N, T >::nnz( const T& threshold_ ) const
2424 {
2425  size_t counter = 0;
2426 
2427  const_iterator it = begin(),
2428  it_end = end();
2429  for( ; it != it_end; ++it)
2430  {
2431  if ( fabs(*it) > threshold_ ) {
2432  ++counter;
2433  }
2434  }
2435 
2436  return counter;
2437 }
2438 
2439 template< size_t M, size_t N, typename T >
2440 void
2441 matrix< M, N, T >::threshold( const T& threshold_value_ )
2442 {
2443  iterator it = begin(),
2444  it_end = end();
2445  for( ; it != it_end; ++it)
2446  {
2447  if ( fabs(*it) <= threshold_value_ ) {
2448  *it = static_cast<T> (0);
2449  }
2450  }
2451 }
2452 
2453 
2454 template< size_t M, size_t N, typename T >
2455 template< typename TT >
2456 void
2457 matrix< M, N, T >::quantize_to( matrix< M, N, TT >& quantized_, const T& min_value, const T& max_value ) const
2458 {
2459  long max_tt_range = long(std::numeric_limits< TT >::max());
2460  long min_tt_range = long(std::numeric_limits< TT >::min());
2461  long tt_range = (max_tt_range - min_tt_range);
2462 
2463  T t_range = max_value - min_value;
2464 
2465  typedef matrix< M, N, TT > m_tt_type ;
2466  typedef typename m_tt_type::iterator tt_iterator;
2467  tt_iterator it_quant = quantized_.begin();
2468  const_iterator it = begin(), it_end = end();
2469 
2470  for( ; it != it_end; ++it, ++it_quant )
2471  {
2472  if (std::numeric_limits<TT>::is_signed ) {
2473  *it_quant = TT( std::min( std::max( min_tt_range, long(( *it * tt_range / t_range ) + 0.5)), max_tt_range ));
2474  } else {
2475  *it_quant = TT( std::min( std::max( min_tt_range, long(((*it - min_value) * tt_range / t_range) + 0.5)), max_tt_range ));
2476  }
2477  }
2478 }
2479 
2480 
2481 template< size_t M, size_t N, typename T >
2482 template< typename TT >
2483 void
2484 matrix< M, N, T >::quantize( matrix< M, N, TT >& quantized_, T& min_value, T& max_value ) const
2485 {
2486  min_value = get_min();
2487  max_value = get_max();
2488  quantize_to( quantized_, min_value, max_value );
2489 }
2490 
2491 
2492 
2493 template< size_t M, size_t N, typename T >
2494 template< typename TT >
2495 void
2496 matrix< M, N, T >::dequantize( matrix< M, N, TT >& dequantized_, const TT& min_value, const TT& max_value ) const
2497 {
2498  long max_t_range = long(std::numeric_limits< T >::max());
2499  long min_t_range = long(std::numeric_limits< T >::min());
2500  long t_range = (max_t_range - min_t_range);
2501 
2502  TT tt_range = max_value - min_value;
2503 
2504  typedef matrix< M, N, TT > m_tt_type ;
2505  typedef typename m_tt_type::iterator tt_iterator;
2506  tt_iterator it_dequant = dequantized_.begin();
2507  const_iterator it = begin(), it_end = end();
2508  for( ; it != it_end; ++it, ++it_dequant )
2509  {
2510  if (std::numeric_limits<T>::is_signed ) {
2511  *it_dequant = std::min( std::max( min_value, TT((TT(*it) / t_range) * tt_range)), max_value );
2512  } else {
2513  *it_dequant = std::min( std::max( min_value, TT((((TT(*it) / t_range)) * tt_range ) + min_value)), max_value );
2514  }
2515  }
2516 }
2517 
2518 template< size_t M, size_t N, typename T >
2519 void
2520 matrix< M, N, T >::columnwise_sum( vector< N, T>& summed_columns_ ) const
2521 {
2522 
2523  for ( size_t n = 0; n < N; ++n )
2524  {
2525  T value = 0;
2526  for ( size_t m = 0; m < M; ++m )
2527  {
2528  value += at( m, n );
2529  }
2530  summed_columns_.at( n ) = value;
2531  }
2532 }
2533 
2534 template< size_t M, size_t N, typename T >
2535 double
2536 matrix< M, N, T >::sum_elements( ) const
2537 {
2538  double sum = 0.0;
2539 
2540  const_iterator it = begin(), it_end = end();
2541  for( ; it != it_end; ++it )
2542  {
2543  sum += *it;
2544  }
2545 
2546  return sum;
2547 }
2548 
2549 template< size_t M, size_t N, typename T >
2550 template< size_t R>
2551 typename enable_if< R == M && R == N>::type*
2552 matrix< M, N, T >::diag( const vector< R, T >& diag_values_ )
2553 {
2554  zero();
2555  for( size_t r = 0; r < R; ++r )
2556  {
2557  at(r, r) = static_cast< T >( diag_values_.at(r) );
2558  }
2559  return 0;
2560 }
2561 
2562 
2563 template< size_t M, size_t N, typename T >
2564 void
2565 matrix< M, N, T >::set_dct()
2566 {
2567  const double num_rows = M;
2568  double fill_value = 0.0f;
2569  for( size_t row = 0; row < M; ++row )
2570  {
2571  // to receive orthonormality
2572  const double weight = ( row == 0.0 ) ? std::sqrt(1/num_rows) :
2573  std::sqrt(2/num_rows);
2574  for( size_t col = 0; col < N; ++col )
2575  {
2576  fill_value = (2 * col + 1) * row * M_PI / (2*M);
2577  fill_value = std::cos( fill_value );
2578  fill_value *= weight;
2579  at( row, col ) = static_cast< T >( fill_value ) ;
2580  }
2581  }
2582 }
2583 
2584 
2585 template< size_t M, size_t N, typename T >
2586 void
2587 matrix< M, N, T >::set_random( int seed )
2588 {
2589  if ( seed >= 0 )
2590  srand( seed );
2591 
2592  double fillValue = 0.0f;
2593  for( size_t row = 0; row < M; ++row )
2594  {
2595  for( size_t col = 0; col < N; ++col )
2596  {
2597  fillValue = rand();
2598  fillValue /= RAND_MAX;
2599  at( row, col ) = -1.0 + 2.0 * static_cast< double >( fillValue ) ;
2600  }
2601  }
2602 }
2603 
2604 template< size_t M, size_t N, typename T >
2605 void
2606 matrix< M, N, T >::write_to_raw( const std::string& dir_, const std::string& filename_ ) const
2607 {
2608  int dir_length = dir_.size() -1;
2609  int last_separator = dir_.find_last_of( "/");
2610  std::string path = dir_;
2611  if (last_separator < dir_length ) {
2612  path.append( "/" );
2613  }
2614  path.append( filename_ );
2615  //check for format
2616  if( filename_.find( "raw", filename_.size() -3) == std::string::npos )
2617  {
2618  path.append( ".");
2619  path.append( "raw" );
2620  }
2621  std::string path_raw = path;
2622 
2623  std::ofstream outfile;
2624  outfile.open( path_raw.c_str() );
2625  if( outfile.is_open() ) {
2626  size_t len_slice = sizeof(T) * M*N;
2627  outfile.write( (char*)&(*this), len_slice );
2628  outfile.close();
2629  } else {
2630  std::cout << "no file open" << std::endl;
2631  }
2632 }
2633 
2634 template< size_t M, size_t N, typename T >
2635 void
2636 matrix< M, N, T >::read_from_raw( const std::string& dir_, const std::string& filename_ )
2637 {
2638  int dir_length = dir_.size() -1;
2639  int last_separator = dir_.find_last_of( "/");
2640  std::string path = dir_;
2641  if (last_separator < dir_length ) {
2642  path.append( "/" );
2643  }
2644  path.append( filename_ );
2645 
2646  size_t max_file_len = 2147483648u - sizeof(T) ;
2647  size_t len_data = sizeof(T) * size();
2648  char* data = new char[ len_data ];
2649  std::ifstream infile;
2650  infile.open( path.c_str(), std::ios::in);
2651 
2652  if( infile.is_open())
2653  {
2654  iterator it = begin(),
2655  it_end = end();
2656 
2657  while ( len_data > 0 )
2658  {
2659  size_t len_read = (len_data % max_file_len ) > 0 ?
2660  len_data % max_file_len : len_data;
2661  len_data -= len_read;
2662  infile.read( data, len_read );
2663 
2664  T* T_ptr = (T*)&(data[0]);
2665  for( ; (it != it_end) && (len_read > 0); ++it, len_read -= sizeof(T) )
2666  {
2667  *it = *T_ptr; ++T_ptr;
2668  }
2669  }
2670 
2671  delete[] data;
2672  infile.close();
2673  } else {
2674  std::cout << "no file open" << std::endl;
2675  }
2676 
2677 }
2678 
2679 
2680 
2681 
2682 template< size_t M, size_t N, typename T >
2683 void
2684 matrix< M, N, T >::write_csv_file( const std::string& dir_, const std::string& filename_ ) const
2685 {
2686  int dir_length = dir_.size() -1;
2687  int last_separator = dir_.find_last_of( "/");
2688  std::string path = dir_;
2689  if (last_separator < dir_length ) {
2690  path.append( "/" );
2691  }
2692  path.append( filename_ );
2693  //check for format
2694  int suffix_pos = filename_.find( "csv", filename_.size() -3);
2695  if( suffix_pos == (-1)) {
2696  path.append( ".");
2697  path.append( "csv" );
2698  }
2699 
2700  std::ofstream outfile;
2701  outfile.open( path.c_str() );
2702  if( outfile.is_open() ) {
2703  outfile << *this << std::endl;
2704  outfile.close();
2705  } else {
2706  std::cout << "no file open" << std::endl;
2707  }
2708 
2709 }
2710 
2711 template< size_t M, size_t N, typename T >
2712 void
2713 matrix< M, N, T >::read_csv_file( const std::string& dir_, const std::string& filename_ )
2714 {
2715  int dir_length = dir_.size() -1;
2716  int last_separator = dir_.find_last_of( "/");
2717  std::string path = dir_;
2718  if (last_separator < dir_length ) {
2719  path.append( "/" );
2720  }
2721  path.append( filename_ );
2722  //check for format
2723  int suffix_pos = filename_.find( "csv", filename_.size() -3);
2724  if( suffix_pos == (-1)) {
2725  path.append( ".");
2726  path.append( "csv" );
2727  }
2728 
2729  std::ifstream infile;
2730  infile.open( path.c_str(), std::ios::in);
2731  if( infile.is_open() ) {
2732  //TODO: not yet implemented
2733  //infile >> *this >> std::endl;
2734  infile.close();
2735  } else {
2736  std::cout << "no file open" << std::endl;
2737  }
2738 
2739 }
2740 
2741 
2742 template< size_t M, size_t N, typename T >
2743 void
2744 matrix< M, N, T >::sum_rows( matrix< M/2, N, T>& other ) const
2745 {
2746  typedef vector< N, T > row_type;
2747 
2748  row_type* row0 = new row_type;
2749  row_type* row1 = new row_type;
2750 
2751  other.zero();
2752 
2753  for ( size_t row = 0; row < M; ++row )
2754  {
2755  get_row( row++, *row0 );
2756  if ( row < M )
2757  {
2758  get_row( row, *row1 );
2759  *row0 += *row1;
2760  other.set_row( row/2 , *row0 );
2761  }
2762  }
2763 
2764  delete row0;
2765  delete row1;
2766 }
2767 
2768 template< size_t M, size_t N, typename T >
2769 void
2770 matrix< M, N, T >::sum_columns( matrix< M, N/2, T>& other ) const
2771 {
2772  typedef vector< M, T > col_type;
2773 
2774  col_type* col0 = new col_type;
2775  col_type* col1 = new col_type;
2776 
2777  other.zero();
2778 
2779  for ( size_t col = 0; col< N; ++col )
2780  {
2781  get_column( col++, *col0 );
2782  if ( col < N )
2783  {
2784  get_column( col, *col1 );
2785  *col0 += *col1;
2786  other.set_column( col/2, *col0 );
2787  }
2788  }
2789 
2790  delete col0;
2791  delete col1;
2792 }
2793 
2794 
2795 
2796 } // namespace vmml
2797 
2798 #endif
vmml::enable_if
Definition: enable_if.hpp:45
vmml::matrix::matrix
matrix()
Construct a zero-initialized matrix.
Definition: matrix.hpp:823
vmml
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