core/vnl/vnl_matrix_fixed.h
Go to the documentation of this file.
00001 // This is core/vnl/vnl_matrix_fixed.h
00002 #ifndef vnl_matrix_fixed_h_
00003 #define vnl_matrix_fixed_h_
00004 #ifdef VCL_NEEDS_PRAGMA_INTERFACE
00005 #pragma interface
00006 #endif
00007 //:
00008 // \file
00009 // \brief fixed size matrix
00010 //
00011 // \author Andrew W. Fitzgibbon, Oxford RRG
00012 // \date   04 Aug 96
00013 //
00014 // \verbatim
00015 //  Modifications
00016 //   Peter Vanroose, 23 Nov 1996:  added explicit copy constructor
00017 //   LSB (Manchester) 15/03/2001:  added Binary I/O and tidied up the documentation
00018 //   Feb.2002 - Peter Vanroose - brief doxygen comment placed on single line
00019 //   Oct.2002 - Amitha Perera  - separated vnl_matrix and vnl_matrix_fixed,
00020 //                               removed necessity for vnl_matrix_fixed_ref
00021 //   Oct.2002 - Peter Vanroose - added inplace_transpose() method
00022 //   Jul.2003 - Paul Smyth     - fixed end() bug, made op*=() more general
00023 //   Mar.2009 - Peter Vanroose - added arg_min() and arg_max()
00024 //   Oct.2010 - Peter Vanroose - mutators and filling methods now return *this
00025 //   Jan.2011 - Peter Vanroose - added methods set_diagonal() & get_diagonal()
00026 // \endverbatim
00027 //-----------------------------------------------------------------------------
00028 
00029 #include <vcl_cstring.h> // memcpy()
00030 #include <vcl_cassert.h>
00031 #include <vcl_iosfwd.h>
00032 
00033 #include "vnl_matrix.h"
00034 #include "vnl_matrix_ref.h"
00035 #include <vnl/vnl_vector.h>
00036 #include <vnl/vnl_vector_fixed.h> // needed for e.g. vnl_matrix_fixed_mat_vec_mult()
00037 #include <vnl/vnl_c_vector.h>
00038 #include <vnl/vnl_config.h> // for VNL_CONFIG_CHECK_BOUNDS
00039 
00040 export template <class T, unsigned int num_rows, unsigned int num_cols> class vnl_matrix_fixed;
00041 
00042 // This mess is for a MSVC6 workaround.
00043 //
00044 // The problem: the matrix-matrix operator* should be written as a
00045 // non-member function since vxl (currently) forbids the use of member
00046 // templates. However, when declared as
00047 // \code
00048 //     template <class T, unsigned m, unsigned n, unsigned o>
00049 //     matrix<T,m,o> operator*( matrix<T,m,n>, matrix<T,n,o> );
00050 // \endcode
00051 // MSVC6 does not find it. A solution is to declare it as a member
00052 // template. However, the obvious
00053 // \code
00054 //     template <unsigned o>
00055 //     matrix<T,num_rows,o> operator*( matrix<T,num_cols,o> );
00056 // \endcode
00057 // causes an internal compiler error. It turns out that if the new
00058 // template parameter "o" comes _first_, then all is okay. Now, we
00059 // can't change the signature of vnl_matrix_fixed to <unsigned num_cols,
00060 // unsigned num_rows, type>, so we use a "hidden" helper matrix. Except
00061 // that user defined conversion operators and conversion constructors
00062 // are not called for templated functions. So we have to use a helper
00063 // base class. The base class is empty, which means that there is no
00064 // loss in space or time efficiency. Finally, we have:
00065 // \code
00066 //   template <unsigned num_cols, unsigned num_rows, class T>
00067 //   class fake_base { };
00068 //
00069 //   template <class T, unsigned num_rows, unsigned num_cols>
00070 //   class matrix : public fake_base<num_cols,num_rows,T>
00071 //   {
00072 //      template <unsigned o>
00073 //      matrix<T,num_rows,o>  operator*( fake_base<o,num_cols,T> );
00074 //   };
00075 // \endcode
00076 // Notice how "o" is first in the list of template parameters. Since
00077 // base class conversions _are_ performed during template matching,
00078 // matrix<T,m,n> is matched as fake_base<n,m,T>, and all is good. For
00079 // some values of good.
00080 //
00081 // Of course, all this trickery is pre-processed away for conforming
00082 // compilers.
00083 //
00084 template <class T, unsigned int num_rows, unsigned int num_cols>
00085 class vnl_matrix_fixed;
00086 template <class T, unsigned M, unsigned N>
00087 inline
00088 vnl_vector_fixed<T, M> vnl_matrix_fixed_mat_vec_mult(const vnl_matrix_fixed<T, M, N>& a, const vnl_vector_fixed<T, N>& b);
00089 template <class T, unsigned M, unsigned N, unsigned O>
00090 inline
00091 vnl_matrix_fixed<T, M, O> vnl_matrix_fixed_mat_mat_mult(const vnl_matrix_fixed<T, M, N>& a, const vnl_matrix_fixed<T, N, O>& b);
00092 #ifdef VCL_VC_6
00093 template <unsigned num_cols, unsigned num_rows, class T>
00094 class vnl_matrix_fixed_fake_base
00095 {
00096 };
00097 
00098 #define VNL_MATRIX_FIXED_VCL60_WORKAROUND : public vnl_matrix_fixed_fake_base<num_cols,num_rows,T>
00099 #else
00100 #define VNL_MATRIX_FIXED_VCL60_WORKAROUND /* no workaround. Phew. */
00101 #endif
00102 
00103 //: Fixed size, stack-stored, space-efficient matrix.
00104 // vnl_matrix_fixed is a fixed-length, stack storage vector. It has
00105 // the same storage size as a C-style array. It is not related via
00106 // inheritance to vnl_matrix. However, it can be converted cheaply to
00107 // a vnl_matrix_ref.
00108 //
00109 // Read the overview documentation of vnl_vector_fixed.
00110 // The text there applies here.
00111 template <class T, unsigned int num_rows, unsigned int num_cols>
00112 class vnl_matrix_fixed  VNL_MATRIX_FIXED_VCL60_WORKAROUND
00113 {
00114   T data_[num_rows][num_cols]; // Local storage
00115 
00116  public:
00117   typedef vnl_matrix_fixed<T,num_rows,num_cols> self;
00118   typedef unsigned int size_type;
00119 
00120   //: Construct an empty num_rows*num_cols matrix
00121   vnl_matrix_fixed() {}
00122 
00123   //: Construct an empty num_rows*num_cols matrix
00124   //
00125   // The sole purpose of this constructor is to match the interface of
00126   // vnl_matrix, so that algorithms can template over the matrix type
00127   // itself.  It is illegal to call this constructor without
00128   // <tt>n==num_rows</tt> and <tt>m==num_cols</tt>.
00129   vnl_matrix_fixed( unsigned n, unsigned m )
00130   {
00131     assert( n == num_rows && m == num_cols );
00132   }
00133 
00134   //: Construct an m*n matrix and fill with value
00135   explicit vnl_matrix_fixed(T value)
00136   {
00137     T* p = data_[0];
00138     unsigned int n = num_rows * num_cols;
00139     while (n--)
00140       *p++ = value;
00141   }
00142 
00143   //: Construct an m*n Matrix and copy data into it row-wise.
00144   explicit vnl_matrix_fixed(const T* datablck)
00145   {
00146     vcl_memcpy(data_[0], datablck, num_rows*num_cols*sizeof(T));
00147   }
00148 
00149   //: Construct an m*n Matrix and copy rhs into it.
00150   //  Abort if rhs is not the same size.
00151   vnl_matrix_fixed(const vnl_matrix_fixed& rhs)
00152   {
00153     vcl_memcpy(data_[0], rhs.data_block(), num_rows*num_cols*sizeof(T));
00154   }
00155 
00156   //: Construct an m*n Matrix and copy rhs into it.
00157   //  Abort if rhs is not the same size.
00158   vnl_matrix_fixed(const vnl_matrix<T>& rhs)
00159   {
00160     assert(rhs.rows() == num_rows && rhs.columns() == num_cols);
00161     vcl_memcpy(data_[0], rhs.data_block(), num_rows*num_cols*sizeof(T));
00162   }
00163 
00164   //  Destruct the m*n matrix.
00165   // An explicit destructor seems to be necessary, at least for gcc 3.0.0,
00166   // to avoid the compiler generating multiple versions of it.
00167   // (This way, a weak symbol is generated; otherwise not.  A bug of gcc 3.0.)
00168   ~vnl_matrix_fixed() {}
00169 
00170   //: Set all elements to value v
00171   // Complexity $O(r.c)$
00172   vnl_matrix_fixed& operator= (T const&v) { fill(v); return *this; }
00173 
00174   //: Copy a vnl_matrix into this.
00175   //  Abort if rhs is not the same size.
00176   vnl_matrix_fixed& operator=(const vnl_matrix<T>& rhs)
00177   {
00178     assert(rhs.rows() == num_rows && rhs.columns() == num_cols);
00179     vcl_memcpy(data_[0], rhs.data_block(), num_rows*num_cols*sizeof(T));
00180     return *this;
00181   }
00182 
00183   //: Copy another vnl_matrix_fixed<T,m,n> into this.
00184   vnl_matrix_fixed& operator=(const vnl_matrix_fixed& rhs)
00185   {
00186     vcl_memcpy(data_[0], rhs.data_block(), num_rows*num_cols*sizeof(T));
00187     return *this;
00188   }
00189 
00190 // Basic 2D-Array functionality-------------------------------------------
00191 
00192   //: Return number of rows
00193   unsigned rows()    const { return num_rows; }
00194 
00195   //: Return number of columns
00196   // A synonym for cols()
00197   unsigned columns()  const { return num_cols; }
00198 
00199   //: Return number of columns
00200   // A synonym for columns()
00201   unsigned cols()    const { return num_cols; }
00202 
00203   //: Return number of elements
00204   // This equals rows() * cols()
00205   unsigned size()    const { return num_rows*num_cols; }
00206 
00207   //: set element
00208   void put (unsigned r, unsigned c, T const& v) { (*this)(r,c) = v; }
00209 
00210   //: set element, and return *this
00211   vnl_matrix_fixed& set (unsigned r, unsigned c, T const& v) { (*this)(r,c) = v; return *this; }
00212 
00213   //: get element
00214   T    get (unsigned r, unsigned c) const { return (*this)(r,c); }
00215 
00216   //: return pointer to given row
00217   // No boundary checking here.
00218   T       * operator[] (unsigned r) { return data_[r]; }
00219 
00220   //: return pointer to given row
00221   // No boundary checking here.
00222   T const * operator[] (unsigned r) const { return data_[r]; }
00223 
00224   //: Access an element for reading or writing
00225   // There are assert style boundary checks - #define NDEBUG to turn them off.
00226   T       & operator() (unsigned r, unsigned c)
00227   {
00228 #if VNL_CONFIG_CHECK_BOUNDS  && (!defined NDEBUG)
00229     assert(r<rows());   // Check the row index is valid
00230     assert(c<cols());   // Check the column index is valid
00231 #endif
00232     return this->data_[r][c];
00233   }
00234 
00235   //: Access an element for reading
00236   // There are assert style boundary checks - #define NDEBUG to turn them off.
00237   T const & operator() (unsigned r, unsigned c) const
00238   {
00239 #if VNL_CONFIG_CHECK_BOUNDS  && (!defined NDEBUG)
00240     assert(r<rows());   // Check the row index is valid
00241     assert(c<cols());   // Check the column index is valid
00242 #endif
00243     return this->data_[r][c];
00244   }
00245 
00246   // ----------------------- Filling and copying -----------------------
00247 
00248   //: Sets all elements of matrix to specified value, and returns "*this".
00249   //  Complexity $O(r.c)$
00250   //  Returning "*this" allows "chaining" two or more operations:
00251   //  e.g., to set a matrix to a column-normalized all-elements-equal matrix, say
00252   //  \code
00253   //     M.fill(1).normalize_columns();
00254   //  \endcode
00255   //  Returning "*this" also allows passing such a matrix as argument
00256   //  to a function f, without having to name the constructed matrix:
00257   //  \code
00258   //     f(vnl_matrix_fixed<double,5,5>(1.0).normalize_columns());
00259   //  \endcode
00260   vnl_matrix_fixed& fill(T);
00261 
00262   //: Sets all diagonal elements of matrix to specified value; returns "*this".
00263   //  Complexity $O(\min(r,c))$
00264   //  Returning "*this" allows "chaining" two or more operations:
00265   //  e.g., to set a 3x3 matrix to [5 0 0][0 10 0][0 0 15], just say
00266   //  \code
00267   //     M.fill_diagonal(5).scale_row(1,2).scale_column(2,3);
00268   //  \endcode
00269   //  Returning "*this" also allows passing a diagonal-filled matrix as argument
00270   //  to a function f, without having to name the constructed matrix:
00271   //  \code
00272   //     f(vnl_matrix_fixed<double,3,3>().fill_diagonal(5));
00273   //  \endcode
00274   vnl_matrix_fixed& fill_diagonal(T);
00275 
00276   //: Sets the diagonal elements of this matrix to the specified list of values.
00277   //  Returning "*this" allows "chaining" two or more operations: see the
00278   //  reasoning (and the examples) in the documentation for method
00279   //  fill_diagonal().
00280   vnl_matrix_fixed& set_diagonal(vnl_vector<T> const&);
00281 
00282   //: Fills (laminates) this matrix with the given data, then returns it.
00283   //  We assume that the argument points to a contiguous rows*cols array, stored rowwise.
00284   //  No bounds checking on the array.
00285   //  Returning "*this" allows "chaining" two or more operations:
00286   //  e.g., to fill a square matrix column-wise, fill it rowwise then transpose:
00287   //  \code
00288   //     M.copy_in(array).inplace_transpose();
00289   //  \endcode
00290   //  Returning "*this" also allows passing a filled-in matrix as argument
00291   //  to a function f, without having to name the constructed matrix:
00292   //  \code
00293   //     f(vnl_matrix_fixed<double,3,3>().copy_in(array));
00294   //  \endcode
00295   vnl_matrix_fixed& copy_in(T const *);
00296 
00297   //: Fills (laminates) this matrix with the given data, then returns it.
00298   //  A synonym for copy_in()
00299   vnl_matrix_fixed& set(T const *d) { return copy_in(d); }
00300 
00301   //: Fills the given array with this matrix.
00302   //  We assume that the argument points to a contiguous rows*cols array, stored rowwise.
00303   //  No bounds checking on the array.
00304   void copy_out(T *) const;
00305 
00306   //: Transposes this matrix efficiently, if it is square, and returns it.
00307   //  Returning "*this" allows "chaining" two or more operations:
00308   //  e.g., to fill a square matrix column-wise, fill it rowwise then transpose:
00309   //  \code
00310   //     M.copy_in(array).inplace_transpose();
00311   //  \endcode
00312   vnl_matrix_fixed& inplace_transpose();
00313 
00314   // ----------------------- Arithmetic --------------------------------
00315   // note that these functions should not pass scalar as a const&.
00316   // Look what would happen to A /= A(0,0).
00317 
00318   //: Add \a s to each element of lhs matrix in situ
00319   vnl_matrix_fixed& operator+= (T s)
00320   {
00321     self::add( data_block(), s, data_block() ); return *this;
00322   }
00323 
00324   //: Subtract \a s from each element of lhs matrix in situ
00325   vnl_matrix_fixed& operator-= (T s)
00326   {
00327     self::sub( data_block(), s, data_block() ); return *this;
00328   }
00329 
00330   //:
00331   vnl_matrix_fixed& operator*= (T s)
00332   {
00333     self::mul( data_block(), s, data_block() ); return *this;
00334   }
00335 
00336   //:
00337   vnl_matrix_fixed& operator/= (T s)
00338   {
00339     self::div( data_block(), s, data_block() ); return *this;
00340   }
00341 
00342   //:
00343   vnl_matrix_fixed& operator+= (vnl_matrix_fixed const& m)
00344   {
00345     self::add( data_block(), m.data_block(), data_block() ); return *this;
00346   }
00347 
00348   //:
00349   vnl_matrix_fixed& operator+= (vnl_matrix<T> const& m)
00350   {
00351     assert( m.rows() == rows() && m.cols() == cols() );
00352     self::add( data_block(), m.data_block(), data_block() ); return *this;
00353   }
00354 
00355   //:
00356   vnl_matrix_fixed& operator-= (vnl_matrix_fixed const& m)
00357   {
00358     self::sub( data_block(), m.data_block(), data_block() ); return *this;
00359   }
00360 
00361   //:
00362   vnl_matrix_fixed& operator-= (vnl_matrix<T> const& m)
00363   {
00364     assert( m.rows() == rows() && m.cols() == cols() );
00365     self::sub( data_block(), m.data_block(), data_block() );
00366     return *this;
00367   }
00368 
00369   //: Negate all elements of matrix
00370   vnl_matrix_fixed operator- () const
00371   {
00372     vnl_matrix_fixed r;
00373     self::sub( T(0), data_block(), r.data_block() );
00374     return r;
00375   }
00376 
00377   //:
00378   vnl_matrix_fixed& operator*= (vnl_matrix_fixed<T,num_cols,num_cols> const& s)
00379   {
00380     vnl_matrix_fixed<T, num_rows, num_cols> out;
00381     for (unsigned i = 0; i < num_rows; ++i)
00382       for (unsigned j = 0; j < num_cols; ++j)
00383       {
00384         T accum = this->data_[i][0] * s(0,j);
00385         for (unsigned k = 1; k < num_cols; ++k)
00386           accum += this->data_[i][k] * s(k,j);
00387         out(i,j) = accum;
00388       }
00389     return *this = out;
00390   }
00391 
00392 #ifdef VCL_VC_6
00393   template <unsigned o>
00394   vnl_matrix_fixed<T,num_rows,o> operator*( vnl_matrix_fixed_fake_base<o,num_cols,T> const& mat ) const
00395   {
00396     vnl_matrix_fixed<T,num_cols,o> const& b = static_cast<vnl_matrix_fixed<T,num_cols,o> const&>(mat);
00397     return vnl_matrix_fixed_mat_mat_mult<T,num_rows,num_cols,o>( *this, b );
00398   }
00399   vnl_vector_fixed<T, num_rows> operator*( vnl_vector_fixed<T, num_cols> const& b) const
00400   {
00401     return vnl_matrix_fixed_mat_vec_mult<T,num_rows,num_cols>(*this,b);
00402   }
00403 #endif
00404 
00405   ////--------------------------- Additions ----------------------------
00406 
00407   //: Make a new matrix by applying function to each element.
00408   vnl_matrix_fixed apply(T (*f)(T)) const;
00409 
00410   //: Make a new matrix by applying function to each element.
00411   vnl_matrix_fixed apply(T (*f)(T const&)) const;
00412 
00413   //: Return transpose
00414   vnl_matrix_fixed<T,num_cols,num_rows> transpose() const;
00415 
00416   //: Return conjugate transpose
00417   vnl_matrix_fixed<T,num_cols,num_rows> conjugate_transpose() const;
00418 
00419   //: Set values of this matrix to those of M, starting at [top,left]
00420   vnl_matrix_fixed& update(vnl_matrix<T> const&, unsigned top=0, unsigned left=0);
00421 
00422   //: Set the elements of the i'th column to v[i]  (No bounds checking)
00423   vnl_matrix_fixed& set_column(unsigned i, T const * v);
00424 
00425   //: Set the elements of the i'th column to value, then return *this.
00426   vnl_matrix_fixed& set_column(unsigned i, T value );
00427 
00428   //: Set j-th column to v, then return *this.
00429   vnl_matrix_fixed& set_column(unsigned j, vnl_vector<T> const& v);
00430 
00431   //: Set j-th column to v, then return *this.
00432   vnl_matrix_fixed& set_column(unsigned j, vnl_vector_fixed<T,num_rows> const& v);
00433 
00434   //: Set columns to those in M, starting at starting_column, then return *this.
00435   vnl_matrix_fixed& set_columns(unsigned starting_column, vnl_matrix<T> const& M);
00436 
00437   //: Set the elements of the i'th row to v[i]  (No bounds checking)
00438   vnl_matrix_fixed& set_row   (unsigned i, T const * v);
00439 
00440   //: Set the elements of the i'th row to value, then return *this.
00441   vnl_matrix_fixed& set_row   (unsigned i, T value );
00442 
00443   //: Set the i-th row, then return *this.
00444   vnl_matrix_fixed& set_row   (unsigned i, vnl_vector<T> const&);
00445 
00446   //: Set the i-th row, then return *this.
00447   vnl_matrix_fixed& set_row   (unsigned i, vnl_vector_fixed<T,num_cols> const&);
00448 
00449   //: Extract a sub-matrix of size r x c, starting at (top,left)
00450   //  Thus it contains elements  [top,top+r-1][left,left+c-1]
00451   vnl_matrix<T> extract (unsigned r,  unsigned c,
00452                          unsigned top=0, unsigned left=0) const;
00453 
00454   //: Extract a sub-matrix starting at (top,left)
00455   //
00456   //  The output is stored in \a sub_matrix, and it should have the
00457   //  required size on entry.  Thus the result will contain elements
00458   //  [top,top+sub_matrix.rows()-1][left,left+sub_matrix.cols()-1]
00459   void extract ( vnl_matrix<T>& sub_matrix,
00460                  unsigned top=0, unsigned left=0) const;
00461 
00462   //: Get a vector equal to the given row
00463   vnl_vector_fixed<T,num_cols> get_row   (unsigned row) const;
00464 
00465   //: Get a vector equal to the given column
00466   vnl_vector_fixed<T,num_rows> get_column(unsigned col) const;
00467 
00468   //: Get n rows beginning at rowstart
00469   vnl_matrix<T> get_n_rows   (unsigned rowstart, unsigned n) const;
00470 
00471   //: Get n columns beginning at colstart
00472   vnl_matrix<T> get_n_columns(unsigned colstart, unsigned n) const;
00473 
00474   //: Return a vector with the content of the (main) diagonal
00475   vnl_vector<T> get_diagonal() const;
00476 
00477   // ==== mutators ====
00478 
00479   //: Sets this matrix to an identity matrix, then returns "*this".
00480   //  Returning "*this" allows e.g. passing an identity matrix as argument to
00481   //  a function f, without having to name the constructed matrix:
00482   //  \code
00483   //     f(vnl_matrix_fixed<double,5,5>().set_identity());
00484   //  \endcode
00485   //  Returning "*this" also allows "chaining" two or more operations:
00486   //  e.g., to set a 3x3 matrix to [3 0 0][0 2 0][0 0 1], one could say
00487   //  \code
00488   //     M.set_identity().scale_row(0,3).scale_column(1,2);
00489   //  \endcode
00490   //  If the matrix is not square, anyhow set main diagonal to 1, the rest to 0.
00491   vnl_matrix_fixed& set_identity();
00492 
00493   //: Reverses the order of rows, and returns "*this".
00494   //  Returning "*this" allows "chaining" two or more operations:
00495   //  e.g., to flip both up-down and left-right, one could just say
00496   //  \code
00497   //     M.flipud().fliplr();
00498   //  \endcode
00499   vnl_matrix_fixed& flipud();
00500 
00501   //: Reverses the order of columns, and returns "*this".
00502   //  Returning "*this" allows "chaining" two or more operations:
00503   //  e.g., to flip both up-down and left-right, one could just say
00504   //  \code
00505   //     M.flipud().fliplr();
00506   //  \endcode
00507   vnl_matrix_fixed& fliplr();
00508 
00509   //: Normalizes each row so it is a unit vector, and returns "*this".
00510   //  Zero rows are not modified
00511   //  Returning "*this" allows "chaining" two or more operations:
00512   //  e.g., to set a matrix to a row-normalized all-elements-equal matrix, say
00513   //  \code
00514   //     M.fill(1).normalize_rows();
00515   //  \endcode
00516   //  Returning "*this" also allows passing such a matrix as argument
00517   //  to a function f, without having to name the constructed matrix:
00518   //  \code
00519   //     f(vnl_matrix_fixed<double,5,5>(1.0).normalize_rows());
00520   //  \endcode
00521   vnl_matrix_fixed& normalize_rows();
00522 
00523   //: Normalizes each column so it is a unit vector, and returns "*this".
00524   //  Zero columns are not modified
00525   //  Returning "*this" allows "chaining" two or more operations:
00526   //  e.g., to set a matrix to a column-normalized all-elements-equal matrix, say
00527   //  \code
00528   //     M.fill(1).normalize_columns();
00529   //  \endcode
00530   //  Returning "*this" also allows passing such a matrix as argument
00531   //  to a function f, without having to name the constructed matrix:
00532   //  \code
00533   //     f(vnl_matrix_fixed<double,5,5>(1.0).normalize_columns());
00534   //  \endcode
00535   vnl_matrix_fixed& normalize_columns();
00536 
00537   //: Scales elements in given row by a factor T, and returns "*this".
00538   //  Returning "*this" allows "chaining" two or more operations:
00539   //  e.g., to set a 3x3 matrix to [3 0 0][0 2 0][0 0 1], one could say
00540   //  \code
00541   //     M.set_identity().scale_row(0,3).scale_column(1,2);
00542   //  \endcode
00543   vnl_matrix_fixed& scale_row   (unsigned row, T value);
00544 
00545   //: Scales elements in given column by a factor T, and returns "*this".
00546   //  Returning "*this" allows "chaining" two or more operations:
00547   //  e.g., to set a 3x3 matrix to [3 0 0][0 2 0][0 0 1], one could say
00548   //  \code
00549   //     M.set_identity().scale_row(0,3).scale_column(1,2);
00550   //  \endcode
00551   vnl_matrix_fixed& scale_column(unsigned col, T value);
00552 
00553   //: Type def for norms.
00554   typedef typename vnl_c_vector<T>::abs_t abs_t;
00555 
00556   //: Return sum of absolute values of elements
00557   abs_t array_one_norm() const { return vnl_c_vector<T>::one_norm(begin(), size()); }
00558 
00559   //: Return square root of sum of squared absolute element values
00560   abs_t array_two_norm() const { return vnl_c_vector<T>::two_norm(begin(), size()); }
00561 
00562   //: Return largest absolute element value
00563   abs_t array_inf_norm() const { return vnl_c_vector<T>::inf_norm(begin(), size()); }
00564 
00565   //: Return sum of absolute values of elements
00566   abs_t absolute_value_sum() const { return array_one_norm(); }
00567 
00568   //: Return largest absolute value
00569   abs_t absolute_value_max() const { return array_inf_norm(); }
00570 
00571   // $ || M ||_1 := \max_j \sum_i | M_{ij} | $
00572   abs_t operator_one_norm() const;
00573 
00574   // $ || M ||_\inf := \max_i \sum_j | M_{ij} | $
00575   abs_t operator_inf_norm() const;
00576 
00577   //: Return Frobenius norm of matrix (sqrt of sum of squares of its elements)
00578   abs_t frobenius_norm() const { return vnl_c_vector<T>::two_norm(begin(), size()); }
00579 
00580   //: Return Frobenius norm of matrix (sqrt of sum of squares of its elements)
00581   abs_t fro_norm() const { return frobenius_norm(); }
00582 
00583   //: Return RMS of all elements
00584   abs_t rms() const { return vnl_c_vector<T>::rms_norm(begin(), size()); }
00585 
00586   //: Return minimum value of elements
00587   T min_value() const { return vnl_c_vector<T>::min_value(begin(), size()); }
00588 
00589   //: Return maximum value of elements
00590   T max_value() const { return vnl_c_vector<T>::max_value(begin(), size()); }
00591 
00592   //: Return location of minimum value of elements
00593   unsigned arg_min() const { return vnl_c_vector<T>::arg_min(begin(), size()); }
00594 
00595   //: Return location of maximum value of elements
00596   unsigned arg_max() const { return vnl_c_vector<T>::arg_max(begin(), size()); }
00597 
00598   //: Return mean of all matrix elements
00599   T mean() const { return vnl_c_vector<T>::mean(begin(), size()); }
00600 
00601   // predicates
00602 
00603   //: Return true iff the size is zero.
00604   bool empty() const { return num_rows==0 && num_cols==0; }
00605 
00606   //:  Return true if all elements equal to identity.
00607   bool is_identity() const;
00608 
00609   //:  Return true if all elements equal to identity, within given tolerance
00610   bool is_identity(double tol) const;
00611 
00612   //: Return true if all elements equal to zero.
00613   bool is_zero() const;
00614 
00615   //: Return true if all elements equal to zero, within given tolerance
00616   bool is_zero(double tol) const;
00617 
00618   //: Return true if finite
00619   bool is_finite() const;
00620 
00621   //: Return true if matrix contains NaNs
00622   bool has_nans() const;
00623 
00624   //: abort if size is not as expected
00625   // This function does or tests nothing if NDEBUG is defined
00626   void assert_size(unsigned nr_rows, unsigned nr_cols) const
00627   {
00628 #ifndef NDEBUG
00629     assert_size_internal(nr_rows, nr_cols);
00630 #endif
00631   }
00632 
00633   //: abort if matrix contains any INFs or NANs.
00634   // This function does or tests nothing if NDEBUG is defined
00635   void assert_finite() const
00636   {
00637 #ifndef NDEBUG
00638     assert_finite_internal();
00639 #endif
00640   }
00641 
00642   ////----------------------- Input/Output ----------------------------
00643 
00644   // : Read a vnl_matrix from an ascii vcl_istream, automatically determining file size if the input matrix has zero size.
00645   bool read_ascii(vcl_istream& s);
00646 
00647   //--------------------------------------------------------------------------------
00648 
00649   //: Access the contiguous block storing the elements in the matrix row-wise. O(1).
00650   // 1d array, row-major order.
00651   T const* data_block () const { return data_[0]; }
00652 
00653   //: Access the contiguous block storing the elements in the matrix row-wise. O(1).
00654   // 1d array, row-major order.
00655   T      * data_block () { return data_[0]; }
00656 
00657 
00658   //----------------------------------------------------------------------
00659   // Conversion to vnl_matrix_ref.
00660 
00661   // The const version of as_ref should return a const vnl_matrix_ref
00662   // so that the vnl_matrix_ref::non_const() cannot be used on
00663   // it. This prevents a const vnl_matrix_fixed from being cast into a
00664   // non-const vnl_matrix reference, giving a slight increase in type safety.
00665 
00666   //: Explicit conversion to a vnl_matrix_ref.
00667   // This is a cheap conversion for those functions that have an interface
00668   // for vnl_matrix but not for vnl_matrix_fixed. There is also a
00669   // conversion operator that should work most of the time.
00670   // \sa vnl_matrix_ref::non_const
00671   vnl_matrix_ref<T> as_ref() { return vnl_matrix_ref<T>( num_rows, num_cols, data_block() ); }
00672 
00673   //: Explicit conversion to a vnl_matrix_ref.
00674   // This is a cheap conversion for those functions that have an interface
00675   // for vnl_matrix but not for vnl_matrix_fixed. There is also a
00676   // conversion operator that should work most of the time.
00677   // \sa vnl_matrix_ref::non_const
00678   const vnl_matrix_ref<T> as_ref() const { return vnl_matrix_ref<T>( num_rows, num_cols, const_cast<T*>(data_block()) ); }
00679 
00680   //: Cheap conversion to vnl_matrix_ref
00681   // Sometimes, such as with templated functions, the compiler cannot
00682   // use this user-defined conversion. For those cases, use the
00683   // explicit as_ref() method instead.
00684   operator const vnl_matrix_ref<T>() const { return vnl_matrix_ref<T>( num_rows, num_cols, const_cast<T*>(data_block()) ); }
00685 
00686   //: Convert to a vnl_matrix.
00687   const vnl_matrix<T> as_matrix() const { return vnl_matrix<T>(const_cast<T*>(data_block()),num_rows,num_cols); }
00688 
00689   //----------------------------------------------------------------------
00690 
00691   typedef T element_type;
00692 
00693   //: Iterators
00694   typedef T       *iterator;
00695   //: Iterator pointing to start of data
00696   iterator       begin() { return data_[0]; }
00697   //: Iterator pointing to element beyond end of data
00698   iterator       end() { return begin() + size(); }
00699 
00700   //: Const iterators
00701   typedef T const *const_iterator;
00702   //: Iterator pointing to start of data
00703   const_iterator begin() const { return data_[0]; }
00704   //: Iterator pointing to element beyond end of data
00705   const_iterator end() const { return begin() + size(); }
00706 
00707   //--------------------------------------------------------------------------------
00708 
00709   //: Return true if *this == rhs
00710   bool operator_eq (vnl_matrix_fixed const & rhs) const
00711   {
00712     return equal( this->data_block(), rhs.data_block() );
00713   }
00714 
00715   //: Equality operator
00716   bool operator==(vnl_matrix_fixed const &that) const { return  this->operator_eq(that); }
00717 
00718   //: Inequality operator
00719   bool operator!=(vnl_matrix_fixed const &that) const { return !this->operator_eq(that); }
00720 
00721   //: Equality operator
00722   bool operator==(vnl_matrix<T> const &that) const { return  this->operator_eq(that); }
00723 
00724   //: Inequality operator
00725   bool operator!=(vnl_matrix<T> const &that) const { return !this->operator_eq(that); }
00726 
00727   //: Print matrix to os in some hopefully sensible format
00728   void print(vcl_ostream& os) const;
00729 
00730 //--------------------------------------------------------------------------------
00731 
00732 
00733   // Helper routines for arithmetic. These routines know the size from
00734   // the template parameters. The vector-vector operations are
00735   // element-wise.
00736 
00737   static void add( const T* a, const T* b, T* r );
00738   static void add( const T* a, T b, T* r );
00739   static void sub( const T* a, const T* b, T* r );
00740   static void sub( const T* a, T b, T* r );
00741   static void sub( T a, const T* b, T* r );
00742   static void mul( const T* a, const T* b, T* r );
00743   static void mul( const T* a, T b, T* r );
00744   static void div( const T* a, const T* b, T* r );
00745   static void div( const T* a, T b, T* r );
00746 
00747   static bool equal( const T* a, const T* b );
00748 
00749  private:
00750   void assert_finite_internal() const;
00751 
00752   void assert_size_internal(unsigned, unsigned) const;
00753 };
00754 
00755 #undef VNL_MATRIX_FIXED_VCL60_WORKAROUND
00756 
00757 
00758 // Make the operators below inline because (1) they are small and
00759 // (2) we then have less explicit instantiation trouble.
00760 
00761 
00762 // --- Matrix-scalar -------------------------------------------------------------
00763 
00764 template <class T, unsigned m, unsigned n>
00765 inline
00766 vnl_matrix_fixed<T,m,n> operator+( const vnl_matrix_fixed<T,m,n>& mat1, const vnl_matrix_fixed<T,m,n>& mat2 )
00767 {
00768   vnl_matrix_fixed<T,m,n> r;
00769   vnl_matrix_fixed<T,m,n>::add( mat1.data_block(), mat2.data_block(), r.data_block() );
00770   return r;
00771 }
00772 
00773 template <class T, unsigned m, unsigned n>
00774 inline
00775 vnl_matrix_fixed<T,m,n> operator+( const vnl_matrix_fixed<T,m,n>& mat, T s )
00776 {
00777   vnl_matrix_fixed<T,m,n> r;
00778   vnl_matrix_fixed<T,m,n>::add( mat.data_block(), s, r.data_block() );
00779   return r;
00780 }
00781 
00782 template <class T, unsigned m, unsigned n>
00783 inline
00784 vnl_matrix_fixed<T,m,n> operator+( const T& s,
00785                                    const vnl_matrix_fixed<T,m,n>& mat )
00786 {
00787   vnl_matrix_fixed<T,m,n> r;
00788   vnl_matrix_fixed<T,m,n>::add( mat.data_block(), s, r.data_block() );
00789   return r;
00790 }
00791 
00792 template <class T, unsigned m, unsigned n>
00793 inline
00794 vnl_matrix_fixed<T,m,n> operator-( const vnl_matrix_fixed<T,m,n>& mat1, const vnl_matrix_fixed<T,m,n>& mat2 )
00795 {
00796   vnl_matrix_fixed<T,m,n> r;
00797   vnl_matrix_fixed<T,m,n>::sub( mat1.data_block(), mat2.data_block(), r.data_block() );
00798   return r;
00799 }
00800 
00801 template <class T, unsigned m, unsigned n>
00802 inline
00803 vnl_matrix_fixed<T,m,n> operator-( const vnl_matrix_fixed<T,m,n>& mat, T s )
00804 {
00805   vnl_matrix_fixed<T,m,n> r;
00806   vnl_matrix_fixed<T,m,n>::sub( mat.data_block(), s, r.data_block() );
00807   return r;
00808 }
00809 
00810 template <class T, unsigned m, unsigned n>
00811 inline
00812 vnl_matrix_fixed<T,m,n> operator-( const T& s,
00813                                    const vnl_matrix_fixed<T,m,n>& mat )
00814 {
00815   vnl_matrix_fixed<T,m,n> r;
00816   vnl_matrix_fixed<T,m,n>::sub( s, mat.data_block(), r.data_block() );
00817   return r;
00818 }
00819 
00820 template <class T, unsigned m, unsigned n>
00821 inline
00822 vnl_matrix_fixed<T,m,n> operator*( const vnl_matrix_fixed<T,m,n>& mat, T s )
00823 {
00824   vnl_matrix_fixed<T,m,n> r;
00825   vnl_matrix_fixed<T,m,n>::mul( mat.data_block(), s, r.data_block() );
00826   return r;
00827 }
00828 
00829 template <class T, unsigned m, unsigned n>
00830 inline
00831 vnl_matrix_fixed<T,m,n> operator*( const T& s,
00832                                    const vnl_matrix_fixed<T,m,n>& mat )
00833 {
00834   vnl_matrix_fixed<T,m,n> r;
00835   vnl_matrix_fixed<T,m,n>::mul( mat.data_block(), s, r.data_block() );
00836   return r;
00837 }
00838 
00839 template <class T, unsigned m, unsigned n>
00840 inline
00841 vnl_matrix_fixed<T,m,n> operator/( const vnl_matrix_fixed<T,m,n>& mat, T s )
00842 {
00843   vnl_matrix_fixed<T,m,n> r;
00844   vnl_matrix_fixed<T,m,n>::div( mat.data_block(), s, r.data_block() );
00845   return r;
00846 }
00847 
00848 
00849 template <class T, unsigned m, unsigned n>
00850 inline
00851 vnl_matrix_fixed<T,m,n> element_product( const vnl_matrix_fixed<T,m,n>& mat1,
00852                                          const vnl_matrix_fixed<T,m,n>& mat2 )
00853 {
00854   vnl_matrix_fixed<T,m,n> r;
00855   vnl_matrix_fixed<T,m,n>::mul( mat1.data_block(), mat2.data_block(), r.data_block() );
00856   return r;
00857 }
00858 
00859 
00860 template <class T, unsigned m, unsigned n>
00861 inline
00862 vnl_matrix_fixed<T,m,n> element_quotient( const vnl_matrix_fixed<T,m,n>& mat1,
00863                                           const vnl_matrix_fixed<T,m,n>& mat2)
00864 {
00865   vnl_matrix_fixed<T,m,n> r;
00866   vnl_matrix_fixed<T,m,n>::div( mat1.data_block(), mat2.data_block(), r.data_block() );
00867   return r;
00868 }
00869 
00870 
00871 // The following two functions are helper functions keep the
00872 // matrix-matrix and matrix-vector multiplication code in one place,
00873 // so that bug fixes, etc, can be localized.
00874 template <class T, unsigned M, unsigned N>
00875 inline
00876 vnl_vector_fixed<T, M>
00877 vnl_matrix_fixed_mat_vec_mult(const vnl_matrix_fixed<T, M, N>& a,
00878                               const vnl_vector_fixed<T, N>& b)
00879 {
00880   vnl_vector_fixed<T, M> out;
00881   for (unsigned i = 0; i < M; ++i)
00882   {
00883     T accum = a(i,0) * b(0);
00884     for (unsigned k = 1; k < N; ++k)
00885       accum += a(i,k) * b(k);
00886     out(i) = accum;
00887   }
00888   return out;
00889 }
00890 
00891 template <class T, unsigned M, unsigned N>
00892 inline
00893 vnl_vector_fixed<T, N>
00894 vnl_matrix_fixed_vec_mat_mult(const vnl_vector_fixed<T, M>& a,
00895                               const vnl_matrix_fixed<T, M, N>& b)
00896 {
00897   vnl_vector_fixed<T, N> out;
00898   for (unsigned i = 0; i < N; ++i)
00899   {
00900     T accum = a(0) * b(0,i);
00901     for (unsigned k = 1; k < M; ++k)
00902       accum += a(k) * b(k,i);
00903     out(i) = accum;
00904   }
00905   return out;
00906 }
00907 
00908 // see comment above
00909 template <class T, unsigned M, unsigned N, unsigned O>
00910 inline
00911 vnl_matrix_fixed<T, M, O>
00912 vnl_matrix_fixed_mat_mat_mult(const vnl_matrix_fixed<T, M, N>& a,
00913                               const vnl_matrix_fixed<T, N, O>& b)
00914 {
00915   vnl_matrix_fixed<T, M, O> out;
00916   for (unsigned i = 0; i < M; ++i)
00917     for (unsigned j = 0; j < O; ++j)
00918     {
00919       T accum = a(i,0) * b(0,j);
00920       for (unsigned k = 1; k < N; ++k)
00921         accum += a(i,k) * b(k,j);
00922       out(i,j) = accum;
00923     }
00924   return out;
00925 }
00926 
00927 #ifndef VCL_VC_6
00928 // The version for correct compilers
00929 
00930 //: Multiply  conformant vnl_matrix_fixed (M x N) and vector_fixed (N)
00931 // \relatesalso vnl_vector_fixed
00932 // \relatesalso vnl_matrix_fixed
00933 template <class T, unsigned M, unsigned N>
00934 inline
00935 vnl_vector_fixed<T, M> operator*(const vnl_matrix_fixed<T, M, N>& a, const vnl_vector_fixed<T, N>& b)
00936 {
00937   return vnl_matrix_fixed_mat_vec_mult(a, b);
00938 }
00939 
00940 //: Multiply  conformant vector_fixed (M) and vnl_matrix_fixed (M x N)
00941 // \relatesalso vnl_vector_fixed
00942 // \relatesalso vnl_matrix_fixed
00943 template <class T, unsigned M, unsigned N>
00944 inline
00945 vnl_vector_fixed<T, N> operator*(const vnl_vector_fixed<T, M>& a, const vnl_matrix_fixed<T, M, N>& b)
00946 {
00947   return vnl_matrix_fixed_vec_mat_mult(a, b);
00948 }
00949 
00950 //: Multiply two conformant vnl_matrix_fixed (M x N) times (N x O)
00951 // \relatesalso vnl_matrix_fixed
00952 template <class T, unsigned M, unsigned N, unsigned O>
00953 inline
00954 vnl_matrix_fixed<T, M, O> operator*(const vnl_matrix_fixed<T, M, N>& a, const vnl_matrix_fixed<T, N, O>& b)
00955 {
00956   return vnl_matrix_fixed_mat_mat_mult(a, b);
00957 }
00958 #endif // VCL_VC_6
00959 
00960 
00961 // These overloads for the common case of mixing a fixed with a
00962 // non-fixed. Because the operator* are templated, the fixed will not
00963 // be automatically converted to a non-fixed-ref. These do it for you.
00964 
00965 template <class T, unsigned m, unsigned n>
00966 inline vnl_matrix<T> operator+( const vnl_matrix_fixed<T,m,n>& a, const vnl_matrix<T>& b )
00967 {
00968   return a.as_ref() + b;
00969 }
00970 
00971 template <class T, unsigned m, unsigned n>
00972 inline vnl_matrix<T> operator+( const vnl_matrix<T>& a, const vnl_matrix_fixed<T,m,n>& b )
00973 {
00974   return a + b.as_ref();
00975 }
00976 
00977 template <class T, unsigned m, unsigned n>
00978 inline vnl_matrix<T> operator-( const vnl_matrix_fixed<T,m,n>& a, const vnl_matrix<T>& b )
00979 {
00980   return a.as_ref() - b;
00981 }
00982 
00983 template <class T, unsigned m, unsigned n>
00984 inline vnl_matrix<T> operator-( const vnl_matrix<T>& a, const vnl_matrix_fixed<T,m,n>& b )
00985 {
00986   return a - b.as_ref();
00987 }
00988 
00989 template <class T, unsigned m, unsigned n>
00990 inline vnl_matrix<T> operator*( const vnl_matrix_fixed<T,m,n>& a, const vnl_matrix<T>& b )
00991 {
00992   return a.as_ref() * b;
00993 }
00994 
00995 template <class T, unsigned m, unsigned n>
00996 inline vnl_matrix<T> operator*( const vnl_matrix<T>& a, const vnl_matrix_fixed<T,m,n>& b )
00997 {
00998   return a * b.as_ref();
00999 }
01000 
01001 template <class T, unsigned m, unsigned n>
01002 inline vnl_vector<T> operator*( const vnl_matrix_fixed<T,m,n>& a, const vnl_vector<T>& b )
01003 {
01004   return a.as_ref() * b;
01005 }
01006 
01007 template <class T, unsigned n>
01008 inline vnl_vector<T> operator*( const vnl_matrix<T>& a, const vnl_vector_fixed<T,n>& b )
01009 {
01010   return a * b.as_ref();
01011 }
01012 
01013 
01014 // --- I/O operations ------------------------------------------------------------
01015 
01016 template <class T, unsigned m, unsigned n>
01017 inline
01018 vcl_ostream& operator<< (vcl_ostream& os, vnl_matrix_fixed<T,m,n> const& mat)
01019 {
01020   mat.print(os);
01021   return os;
01022 }
01023 
01024 template <class T, unsigned m, unsigned n>
01025 inline
01026 vcl_istream& operator>> (vcl_istream& is, vnl_matrix_fixed<T,m,n>& mat)
01027 {
01028   mat.read_ascii(is);
01029   return is;
01030 }
01031 
01032 // More workarounds for Visual C++ 6.0. The problem is that VC6 cannot
01033 // automatically determine the m of the second parameter, for some
01034 // reason. Also, VC6 can't figure out that vector_fixed::SIZE is a
01035 // compile time constant when used in the return parameter. So, we
01036 // have to introduce a helper class to do it.
01037 //
01038 #if defined(VCL_VC_6) && !defined(__GCCXML__)
01039 
01040 template<class T, unsigned m, class FixedVector>
01041 struct outer_product_fixed_type_helper
01042 {
01043   typedef vnl_matrix_fixed<T,m,FixedVector::SIZE> result_matrix;
01044 };
01045 
01046 template<class V1, class V2, class RM>
01047 struct outer_product_fixed_calc_helper
01048 {
01049   static RM calc( V1 const& a, V2 const& b );
01050 };
01051 
01052 template <class T, unsigned m, class SecondFixedVector>
01053 outer_product_fixed_type_helper<T,m,SecondFixedVector>::result_matrix
01054 outer_product(vnl_vector_fixed<T,m> const& a, SecondFixedVector const& b)
01055 {
01056   typedef vnl_vector_fixed<T,m> VecA;
01057   typedef vnl_vector_fixed<T,SecondFixedVector::SIZE> VecB;
01058   typedef outer_product_fixed_type_helper<T,m,SecondFixedVector>::result_matrix ResultMat;
01059   return outer_product_fixed_calc_helper<VecA,VecB,ResultMat>::calc(a,b);
01060 }
01061 
01062 #else // no need for VC6 workaround for outer_product
01063 
01064 //:
01065 // \relatesalso vnl_vector_fixed
01066 template <class T, unsigned m, unsigned n>
01067 vnl_matrix_fixed<T,m,n> outer_product(vnl_vector_fixed<T,m> const& a, vnl_vector_fixed<T,n> const& b);
01068 
01069 #endif // VC6 workaround for outer_product
01070 
01071 #define VNL_MATRIX_FIXED_INSTANTIATE(T, M, N) \
01072 extern "please include vnl/vnl_matrix_fixed.txx instead"
01073 
01074 #endif // vnl_matrix_fixed_h_