core/vnl/vnl_matrix_fixed_ref.h
Go to the documentation of this file.
00001 // This is core/vnl/vnl_matrix_fixed_ref.h
00002 #ifndef vnl_matrix_fixed_ref_h_
00003 #define vnl_matrix_fixed_ref_h_
00004 #ifdef VCL_NEEDS_PRAGMA_INTERFACE
00005 #pragma interface
00006 #endif
00007 //:
00008 // \file
00009 // \brief Fixed size stack-stored vnl_matrix
00010 //
00011 // vnl_matrix_fixed_ref is a fixed-size vnl_matrix for which the data space
00012 // has been supplied externally.  This is useful for two main tasks:
00013 //
00014 // (a) Treating some row-based "C" matrix as a vnl_matrix in order to
00015 // perform vnl_matrix operations on it.
00016 //
00017 // (b) Declaring a vnl_matrix that uses entirely stack-based storage for the
00018 // matrix.
00019 //
00020 // The big warning is that returning a vnl_matrix_fixed_ref pointer will free
00021 // non-heap memory if deleted through a vnl_matrix pointer.  This should be
00022 // very difficult though, as vnl_matrix_fixed_ref objects may not be constructed
00023 // using operator new.  This in turn is plausible as the point is to avoid
00024 // such calls.
00025 //
00026 // \author Andrew W. Fitzgibbon, Oxford RRG
00027 // \date   04 Aug 1996
00028 //
00029 // Additional comments on the vnl_matrix_fixed_ref and vnl_vector_fixed_ref
00030 // classes, extracted from an email conversation between Paul P. Smyth,
00031 // Vicon Motion Systems Ltd., from May 02, 2001, and Amitha Perera
00032 // (who answers the following on Monday, October 07, 2002):
00033 //
00034 // I'm working on separating vnl_vector and vnl_vector_fixed in the VXL
00035 // tree, as I mailed a while ago to the vxl-maintainers list. I noticed
00036 // that you'd committed a vnl_vector_fixed_ref class which doesn't seem
00037 // to provide any additional functionality over vnl_vector_ref. May I
00038 // remove it, or is there some use for it?
00039 //
00040 // Paul Smyth writes:
00041 // The rationale behind it was that I had some (fast) algorithms for
00042 // matrix/vector operations that made use of compile-time knowledge of the
00043 // vector and matrix sizes.
00044 // This was typically appropriate when trying to interpret a fixed-size
00045 // subvector within a large vector of parameters as e.g. a translation.
00046 //
00047 // As I saw it, the various types of vector possible were: (with their current
00048 // names)
00049 // - pointer to memory, plus compile-time knowledge of vector size ( T*, and enum{size}) = vnl_vector_fixed_ref
00050 // - ownership of memory, plus compile-time size = vnl_vector_fixed
00051 // - pointer to memory, plus run-time only knowledge of size (T* and size()) = vnl_vector_ref
00052 // - ownership of memory, variably sized = vnl_vector
00053 //
00054 // I had a conversation with Andrew Fitzgibbon, where he reckoned that the best
00055 // thing to do with vnl vectors etc. was to create entirely separate types, and
00056 // routines for conversion between them (possibly implicitly), rather that
00057 // trying to establish a class hierarchy, which may add too many burdens in
00058 // terms of object size for small vectors/matrices.
00059 //
00060 // Sorry - I've now found the debate on the maintainers list!
00061 //
00062 // Anyway, I believe that vector_fixed_ref is very necessary, and that you
00063 // should be able to convert from a vector_fixed to a vector_fixed_ref - say
00064 // using an as_ref() member on vector_fixed or standalone function.
00065 // And I believe that for the restructured classes, vector_fixed_ref and
00066 // vector_fixed should not be related by inheritance, as that would place an
00067 // excessive burden on the size of vector_fixed.
00068 //
00069 // ------
00070 // Another issue - do you have a mechanism for dealing with const data safely?
00071 // {
00072 //   template<typename T, int n>
00073 //   vnl_vector_fixed_ref(T* i_Data);
00074 //
00075 //   void MyFunction(const vnl_vector<double> & Input)
00076 //   {
00077 //     // take a reference to the first 3 elements of Input
00078 //     vnl_vector_fixed_ref<double,3> ref(Input.begin());
00079 //     // compiler error - as making vector_fixed_ref from const
00080 // double *
00081 //   }
00082 // }
00083 //
00084 // The options appear to be
00085 // 1) Make a separate class vnl_vector_fixed_ref_const
00086 // 2) Make vnl_vector_fixed_ref so it can be instantiated with
00087 // vnl_vector_fixed_ref<double,n> AND vnl_vector_fixed_ref<const double,n>, and
00088 // gives appropriate behaviour - would probably require a to_const function
00089 // which generates vnl_vector_fixed_ref<const T,n> from
00090 // vnl_vector_fixed_ref<T,n>
00091 //
00092 // ------
00093 // Another note is that a number of routines that use vector_fixed currently
00094 // (e.g. cross_3d) should really use vector_fixed_ref as an input, because they
00095 // should be able to operate on fixed vector references as well as fixed
00096 // vectors.
00097 //
00098 // While I'm at it, has it been decided that the vnl_vector and vnl_vector_ref
00099 // classes are to remain unchanged? Because having vnl_vector as the base, and
00100 // vnl_vector_ref derived from it is a real pain in the backside. A vector
00101 // which may or may not own its own memory is a more general type than one
00102 // which does own its own memory, and having vnl_vector as the base means that
00103 // all sorts of nastinesses can happen. Simply, a vector_ref Is-not a type of
00104 // vector.
00105 // If anything, it should be the other way round.
00106 //
00107 // void DoAssign(vnl_vector<double> & RefToMemoryIDontOwn, const vnl_vector<double> & NewContents)
00108 // {
00109 //   RefToMemoryIDontOwn = NewContents;
00110 // }
00111 //
00112 // void DeleteTwice()
00113 // {
00114 //   vnl_vector<double> vec1(3, 0); // size 3 - news 3*double
00115 //   vnl_vector<double> vec2(4,1); // size 4 news 4 * double
00116 //   vnl_vector_ref<double> ref_to_1(3,vec1.begin()); // copies pointer
00117 //   DoAssign(ref_to_1, vec2); // deletes memory owned by 1, news 4 * double
00118 //   // vec1 now points to deleted memory, and will crash when goes out of scope
00119 // }
00120 //
00121 // Maybe that issue isn't on your agenda - but it's a bit of a disaster. I know
00122 // that fixing this might break some code.
00123 //
00124 // ---------
00125 // Sorry for rolling all these things into one - I'd be interested to know what
00126 // you think. But please don't kill my vnl_vector_ref!
00127 //
00128 // Paul.
00129 //
00130 // \verbatim
00131 //  Modifications:
00132 //   27-Nov-1996 Peter Vanroose - added default constructor which allocates matrix storage
00133 //    4-Jul-2003 Paul Smyth - general cleanup and rewrite; interface now as vnl_matrix_fixed
00134 //   15-Aug-2003 Peter Vanroose - removed "duplicate" operator=(vnl_matrix_fixed<T,n> const&)
00135 //    8-Dec-2006 Markus Moll - changed operator>> signature (to const& argument)
00136 //   30-Mar-2009 Peter Vanroose - added arg_min() and arg_max()
00137 //   24-Oct-2010 Peter Vanroose - mutators and filling methods now return *this
00138 //   18-Jan-2011 Peter Vanroose - added methods set_diagonal() & get_diagonal()
00139 // \endverbatim
00140 //
00141 //-----------------------------------------------------------------------------
00142 
00143 #include <vcl_cassert.h>
00144 #include <vcl_iosfwd.h>
00145 #include <vcl_cstring.h> // memcpy()
00146 #include <vnl/vnl_matrix_fixed.h>
00147 #include <vnl/vnl_vector_fixed.h>
00148 #include <vnl/vnl_vector_fixed_ref.h>
00149 #include <vnl/vnl_c_vector.h>
00150 
00151 //: Fixed size stack-stored vnl_matrix
00152 // vnl_matrix_fixed_ref is a fixed-size vnl_matrix for which the data space
00153 // has been supplied externally.  This is useful for two main tasks:
00154 //
00155 // (a) Treating some row-based "C" matrix as a vnl_matrix in order to
00156 // perform vnl_matrix operations on it.
00157 //
00158 // (b) Declaring a vnl_matrix that uses entirely stack-based storage for the
00159 // matrix.
00160 //
00161 // The big warning is that returning a vnl_matrix_fixed_ref pointer will free
00162 // non-heap memory if deleted through a vnl_matrix pointer.  This should be
00163 // very difficult though, as vnl_matrix_fixed_ref objects may not be constructed
00164 // using operator new.  This in turn is plausible as the point is to avoid
00165 // such calls.
00166 //
00167 
00168 template <class T, unsigned num_rows, unsigned num_cols>
00169 class vnl_matrix_fixed_ref_const
00170 {
00171  protected:
00172   const T* data_;
00173  public:
00174   vnl_matrix_fixed_ref_const(const vnl_matrix_fixed<T,num_rows,num_cols>& rhs)
00175     : data_(rhs.data_block())
00176   {
00177   }
00178   explicit vnl_matrix_fixed_ref_const(const T * dataptr)
00179     : data_(dataptr)
00180   {
00181   }
00182   vnl_matrix_fixed_ref_const(const vnl_matrix_fixed_ref_const<T,num_rows,num_cols> & rhs)
00183     : data_(rhs.data_)
00184   {
00185   }
00186   //: Get j-th row
00187   vnl_vector_fixed<T,num_rows> get_row(unsigned row_index) const
00188   {
00189     vnl_vector<T> v(num_cols);
00190     for (unsigned int j = 0; j < num_cols; j++)    // For each element in row
00191       v[j] = (*this)(row_index,j);
00192     return v;
00193   }
00194 
00195   //: Get j-th column
00196   vnl_vector_fixed<T,num_cols> get_column(unsigned column_index) const
00197   {
00198     vnl_vector<T> v(num_rows);
00199     for (unsigned int j = 0; j < num_rows; j++)
00200       v[j] = (*this)(j,column_index);
00201     return v;
00202   }
00203 
00204   //: Return a vector with the content of the (main) diagonal
00205   vnl_vector<T> get_diagonal() const;
00206 
00207   const T * data_block() const { return data_; }
00208 
00209   //: Const iterators
00210   typedef T const *const_iterator;
00211   //: Iterator pointing to start of data
00212   const_iterator begin() const { return data_; }
00213   //: Iterator pointing to element beyond end of data
00214   const_iterator end() const { return begin() + this->size(); }
00215 
00216   //: Type defs for iterators
00217   typedef const T element_type;
00218   //: Type defs for iterators
00219   typedef const T       *iterator;
00220 
00221   T const & operator() (unsigned r, unsigned c) const
00222   {
00223 #if VNL_CONFIG_CHECK_BOUNDS  && (!defined NDEBUG)
00224     assert(r<num_rows);   // Check the row index is valid
00225     assert(c<num_cols);   // Check the column index is valid
00226 #endif
00227     return *(data_ + num_cols * r + c);
00228   }
00229 
00230   //: return pointer to given row
00231   // No boundary checking here.
00232   T const * operator[] (unsigned r) const { return data_ + num_cols * r; }
00233 
00234   //: Return number of rows
00235   unsigned rows()    const { return num_rows; }
00236 
00237   //: Return number of columns
00238   // A synonym for cols()
00239   unsigned columns()  const { return num_cols; }
00240 
00241   //: Return number of columns
00242   // A synonym for columns()
00243   unsigned cols()    const { return num_cols; }
00244 
00245   //: Return number of elements
00246   // This equals rows() * cols()
00247   unsigned size()    const { return num_rows*num_cols; }
00248 
00249   //: Print matrix to os in some hopefully sensible format
00250   void print(vcl_ostream& os) const;
00251 
00252   void copy_out(T *) const;
00253 
00254   ////--------------------------- Additions ----------------------------
00255 
00256   //: Make a new matrix by applying function to each element.
00257   vnl_matrix_fixed<T,num_rows,num_cols> apply(T (*f)(T)) const;
00258 
00259   //: Make a new matrix by applying function to each element.
00260   vnl_matrix_fixed<T,num_rows,num_cols> apply(T (*f)(T const&)) const;
00261 
00262   //: Return transpose
00263   vnl_matrix_fixed<T,num_cols,num_rows> transpose () const;
00264 
00265   //: Return conjugate transpose
00266   vnl_matrix_fixed<T,num_cols,num_rows> conjugate_transpose () const;
00267 
00268   //: Extract a sub-matrix of size rows x cols, starting at (top,left)
00269   //  Thus it contains elements  [top,top+rows-1][left,left+cols-1]
00270   vnl_matrix<T> extract (unsigned rowz,  unsigned colz,
00271                          unsigned top=0, unsigned left=0) const;
00272 
00273   //: Get n rows beginning at rowstart
00274   vnl_matrix<T> get_n_rows   (unsigned rowstart, unsigned n) const;
00275 
00276   //: Get n columns beginning at colstart
00277   vnl_matrix<T> get_n_columns(unsigned colstart, unsigned n) const;
00278 
00279   //: Type def for norms.
00280   typedef typename vnl_c_vector<T>::abs_t abs_t;
00281 
00282   //: Return sum of absolute values of elements
00283   abs_t array_one_norm() const { return vnl_c_vector<T>::one_norm(begin(), size()); }
00284 
00285   //: Return square root of sum of squared absolute element values
00286   abs_t array_two_norm() const { return vnl_c_vector<T>::two_norm(begin(), size()); }
00287 
00288   //: Return largest absolute element value
00289   abs_t array_inf_norm() const { return vnl_c_vector<T>::inf_norm(begin(), size()); }
00290 
00291   //: Return sum of absolute values of elements
00292   abs_t absolute_value_sum() const { return array_one_norm(); }
00293 
00294   //: Return largest absolute value
00295   abs_t absolute_value_max() const { return array_inf_norm(); }
00296 
00297   // $ || M ||_1 := \max_j \sum_i | M_{ij} | $
00298   abs_t operator_one_norm() const;
00299 
00300   // $ || M ||_\inf := \max_i \sum_j | M_{ij} | $
00301   abs_t operator_inf_norm() const;
00302 
00303   //: Return Frobenius norm of matrix (sqrt of sum of squares of its elements)
00304   abs_t frobenius_norm() const { return vnl_c_vector<T>::two_norm(begin(), size()); }
00305 
00306   //: Return Frobenius norm of matrix (sqrt of sum of squares of its elements)
00307   abs_t fro_norm() const { return frobenius_norm(); }
00308 
00309   //: Return RMS of all elements
00310   abs_t rms() const { return vnl_c_vector<T>::rms_norm(begin(), size()); }
00311 
00312   //: Return minimum value of elements
00313   T min_value() const { return vnl_c_vector<T>::min_value(begin(), size()); }
00314 
00315   //: Return maximum value of elements
00316   T max_value() const { return vnl_c_vector<T>::max_value(begin(), size()); }
00317 
00318   //: Return location of minimum value of elements
00319   unsigned arg_min() const { return vnl_c_vector<T>::arg_min(begin(), size()); }
00320 
00321   //: Return location of maximum value of elements
00322   unsigned arg_max() const { return vnl_c_vector<T>::arg_max(begin(), size()); }
00323 
00324   //: Return mean of all matrix elements
00325   T mean() const { return vnl_c_vector<T>::mean(begin(), size()); }
00326 
00327   // predicates
00328 
00329   //: Return true iff the size is zero.
00330   bool empty() const { return num_rows==0 && num_cols==0; }
00331 
00332   //:  Return true if all elements equal to identity.
00333   bool is_identity() const;
00334 
00335   //:  Return true if all elements equal to identity, within given tolerance
00336   bool is_identity(double tol) const;
00337 
00338   //: Return true if all elements equal to zero.
00339   bool is_zero() const;
00340 
00341   //: Return true if all elements equal to zero, within given tolerance
00342   bool is_zero(double tol) const;
00343 
00344   //: Return true if finite
00345   bool is_finite() const;
00346 
00347   //: Return true if matrix contains NaNs
00348   bool has_nans() const;
00349 
00350   //: abort if size is not as expected
00351   // This function does or tests nothing if NDEBUG is defined
00352   void assert_size(unsigned rowz, unsigned colz) const
00353   {
00354 #ifndef NDEBUG
00355     assert_size_internal(rowz, colz);
00356 #endif
00357   }
00358   //: abort if matrix contains any INFs or NANs.
00359   // This function does or tests nothing if NDEBUG is defined
00360   void assert_finite() const
00361   {
00362 #ifndef NDEBUG
00363     assert_finite_internal();
00364 #endif
00365   }
00366 
00367   static void add( const T* a, const T* b, T* r ) { vnl_matrix_fixed<T,num_rows,num_cols>::add(a,b,r); }
00368   static void add( const T* a, T b, T* r ) { vnl_matrix_fixed<T,num_rows,num_cols>::add(a,b,r); }
00369   static void sub( const T* a, const T* b, T* r ) { vnl_matrix_fixed<T,num_rows,num_cols>::sub(a,b,r); }
00370   static void sub( const T* a, T b, T* r ) { vnl_matrix_fixed<T,num_rows,num_cols>::sub(a,b,r); }
00371   static void sub( T a, const T* b, T* r ) { vnl_matrix_fixed<T,num_rows,num_cols>::sub(a,b,r); }
00372   static void mul( const T* a, const T* b, T* r ) { vnl_matrix_fixed<T,num_rows,num_cols>::mul(a,b,r); }
00373   static void mul( const T* a, T b, T* r ) { vnl_matrix_fixed<T,num_rows,num_cols>::mul(a,b,r); }
00374   static void div( const T* a, const T* b, T* r ) { vnl_matrix_fixed<T,num_rows,num_cols>::div(a,b,r); }
00375   static void div( const T* a, T b, T* r ) { vnl_matrix_fixed<T,num_rows,num_cols>::div(a,b,r); }
00376 
00377   static bool equal( const T* a, const T* b ) { return vnl_matrix_fixed<T,num_rows,num_cols>::equal(a,b); }
00378 
00379  private:
00380   const vnl_matrix_fixed_ref_const<T,num_rows,num_cols> & operator=(const vnl_matrix_fixed_ref_const<T,num_rows,num_cols>& ) const
00381   {
00382     assert(!"Assignment is illegal for a vnl_matrix_fixed_ref_const");
00383     return *this;
00384   }
00385 
00386   void assert_finite_internal() const;
00387 
00388   void assert_size_internal(unsigned, unsigned) const;
00389 };
00390 
00391 
00392 template <class T, unsigned num_rows, unsigned num_cols>
00393 class vnl_matrix_fixed_ref : public vnl_matrix_fixed_ref_const<T,num_rows,num_cols>
00394 {
00395   typedef vnl_matrix_fixed_ref_const<T,num_rows,num_cols> base;
00396 
00397  public:
00398   // this is the only point where the const_cast happens
00399   // the base class is used to store the pointer, so that conversion is not necessary
00400   T * data_block() const {
00401     return const_cast<T*>(this->data_);
00402   }
00403   vnl_matrix_fixed_ref(vnl_matrix_fixed<T,num_rows,num_cols>& rhs)
00404     : base(rhs.data_block())
00405   {
00406   }
00407   explicit vnl_matrix_fixed_ref(T * dataptr)
00408     : base(dataptr)
00409   {
00410   }
00411 
00412   //: Copy another vnl_matrix_fixed<T,m,n> into this.
00413   vnl_matrix_fixed_ref const & operator=(const vnl_matrix_fixed_ref_const<T,num_rows,num_cols>& rhs) const
00414   {
00415     vcl_memcpy(data_block(), rhs.data_block(), num_rows*num_cols*sizeof(T));
00416     return *this;
00417   }
00418 
00419   // Basic 2D-Array functionality-------------------------------------------
00420 
00421   //: set element
00422   void put (unsigned r, unsigned c, T const& v) { (*this)(r,c) = v; }
00423 
00424   //: get element
00425   T    get (unsigned r, unsigned c) const { return (*this)(r,c); }
00426 
00427   //: return pointer to given row
00428   // No boundary checking here.
00429   T  * operator[] (unsigned r) const { return data_block() + num_cols * r; }
00430 
00431   //: Access an element for reading or writing
00432   // There are assert style boundary checks - #define NDEBUG to turn them off.
00433   T       & operator() (unsigned r, unsigned c) const
00434   {
00435 #if VNL_CONFIG_CHECK_BOUNDS  && (!defined NDEBUG)
00436     assert(r<num_rows);   // Check the row index is valid
00437     assert(c<num_cols);   // Check the column index is valid
00438 #endif
00439     return *(this->data_block() + num_cols * r + c);
00440   }
00441 
00442   // ----------------------- Filling and copying -----------------------
00443 
00444   //: Sets all elements of matrix to specified value, and returns "*this".
00445   //  Complexity $O(r.c)$
00446   //  Returning "*this" allows "chaining" two or more operations:
00447   //  e.g., to set a matrix to a column-normalized all-elements-equal matrix, say
00448   //  \code
00449   //     M.fill(1).normalize_columns();
00450   //  \endcode
00451   //  Returning "*this" also allows passing such a matrix as argument
00452   //  to a function f, without having to name the constructed matrix:
00453   //  \code
00454   //     f(vnl_matrix_fixed_ref_const<double,5,5>(1.0).normalize_columns());
00455   //  \endcode
00456   vnl_matrix_fixed_ref const& fill (T) const;
00457 
00458   //: Sets all diagonal elements of matrix to specified value; returns "*this".
00459   //  Complexity $O(\min(r,c))$
00460   //  Returning "*this" allows "chaining" two or more operations:
00461   //  e.g., to set a 3x3 matrix to [5 0 0][0 10 0][0 0 15], just say
00462   //  \code
00463   //     M.fill_diagonal(5).scale_row(1,2).scale_column(2,3);
00464   //  \endcode
00465   //  Returning "*this" also allows passing a diagonal-filled matrix as argument
00466   //  to a function f, without having to name the constructed matrix:
00467   //  \code
00468   //     f(vnl_matrix_fixed_ref<double,3,3>().fill_diagonal(5));
00469   //  \endcode
00470   vnl_matrix_fixed_ref const& fill_diagonal (T) const;
00471 
00472   //: Sets the diagonal elements of this matrix to the specified list of values.
00473   //  Returning "*this" allows "chaining" two or more operations: see the
00474   //  reasoning (and the examples) in the documentation for method
00475   //  fill_diagonal().
00476   vnl_matrix_fixed_ref const& set_diagonal(vnl_vector<T> const&) const;
00477 
00478   //: Fills (laminates) this matrix with the given data, then returns it.
00479   //  We assume that the argument points to a contiguous rows*cols array, stored rowwise.
00480   //  No bounds checking on the array.
00481   //  Returning "*this" allows "chaining" two or more operations:
00482   //  e.g., to fill a square matrix column-wise, fill it rowwise then transpose:
00483   //  \code
00484   //     M.copy_in(array).inplace_transpose();
00485   //  \endcode
00486   //  Returning "*this" also allows passing a filled-in matrix as argument
00487   //  to a function f, without having to name the constructed matrix:
00488   //  \code
00489   //     f(vnl_matrix_fixed_ref<double,3,3>().copy_in(array));
00490   //  \endcode
00491   vnl_matrix_fixed_ref const& copy_in(T const *) const;
00492 
00493   //: Fills (laminates) this matrix with the given data, then returns it.
00494   //  A synonym for copy_in()
00495   vnl_matrix_fixed_ref const& set(T const *d) const { return copy_in(d); }
00496 
00497   //: Fills the given array with this matrix.
00498   //  We assume that the argument points to a contiguous rows*cols array, stored rowwise.
00499   //  No bounds checking on the array
00500 
00501   //: Transposes this matrix efficiently, if it is square, and returns it.
00502   //  Returning "*this" allows "chaining" two or more operations:
00503   //  e.g., to fill a square matrix column-wise, fill it rowwise then transpose:
00504   //  \code
00505   //     M.copy_in(array).inplace_transpose();
00506   //  \endcode
00507   vnl_matrix_fixed_ref const& inplace_transpose() const;
00508 
00509   // ----------------------- Arithmetic --------------------------------
00510   // note that these functions should not pass scalar as a const&.
00511   // Look what would happen to A /= A(0,0).
00512 
00513   //: Add \a s to each element of lhs matrix in situ
00514   vnl_matrix_fixed_ref const& operator+= (T s) const
00515   {
00516     base::add( data_block(), s, data_block() ); return *this;
00517   }
00518 
00519   //: Subtract \a s from each element of lhs matrix in situ
00520   vnl_matrix_fixed_ref const& operator-= (T s) const
00521   {
00522     base::sub( data_block(), s, data_block() ); return *this;
00523   }
00524 
00525   //:
00526   vnl_matrix_fixed_ref const& operator*= (T s) const
00527   {
00528     base::mul( data_block(), s, data_block() ); return *this;
00529   }
00530 
00531   //:
00532   vnl_matrix_fixed_ref const& operator/= (T s) const
00533   {
00534     base::div( data_block(), s, data_block() ); return *this;
00535   }
00536 
00537   //:
00538   vnl_matrix_fixed_ref const & operator+= (vnl_matrix_fixed_ref_const<T,num_rows,num_cols> const& m) const
00539   {
00540     base::add( data_block(), m.data_block(), data_block() ); return *this;
00541   }
00542 
00543   //:
00544   vnl_matrix_fixed_ref const& operator+= (vnl_matrix<T> const& m) const
00545   {
00546     assert( m.rows() == num_rows && m.cols() == num_cols );
00547     base::add( data_block(), m.data_block(), data_block() ); return *this;
00548   }
00549 
00550   //:
00551   vnl_matrix_fixed_ref const& operator-= (vnl_matrix_fixed_ref_const<T,num_rows,num_cols> const& m) const
00552   {
00553     base::sub( data_block(), m.data_block(), data_block() ); return *this;
00554   }
00555 
00556   //:
00557   vnl_matrix_fixed_ref const& operator-= (vnl_matrix<T> const& m) const
00558   {
00559     assert( m.rows() == num_rows && m.cols() == num_cols );
00560     base::sub( data_block(), m.data_block(), data_block() ); return *this;
00561   }
00562 
00563   //: Negate all elements of matrix
00564   vnl_matrix_fixed<T,num_rows,num_cols> operator- () const
00565   {
00566     vnl_matrix_fixed<T,num_rows,num_cols> r;
00567     base::sub( T(0), data_block(), r.data_block() );
00568     return r;
00569   }
00570 
00571   //:
00572   vnl_matrix_fixed_ref const& operator*= (vnl_matrix_fixed_ref_const<T,num_cols,num_cols> const& s) const
00573   {
00574     vnl_matrix_fixed<T, num_rows, num_cols> out;
00575     for (unsigned i = 0; i < num_rows; ++i)
00576       for (unsigned j = 0; j < num_cols; ++j)
00577       {
00578         T accum = this->operator()(i,0) * s(0,j);
00579         for (unsigned k = 1; k < num_cols; ++k)
00580           accum += this->operator()(i,k) * s(k,j);
00581         out(i,j) = accum;
00582       }
00583     *this = out;
00584     return *this;
00585   }
00586 
00587 #ifdef VCL_VC_6
00588   template<unsigned o>
00589   vnl_matrix_fixed<T,num_rows,o> operator*( vnl_matrix_fixed_fake_base<o,num_cols,T> const& mat ) const
00590   {
00591     vnl_matrix_fixed<T,num_cols,o> const& b = static_cast<vnl_matrix_fixed<T,num_cols,o> const&>(mat);
00592     return vnl_matrix_fixed_mat_mat_mult<T,num_rows,num_cols,o>( *this, b );
00593   }
00594   vnl_vector_fixed<T, num_rows> operator*( vnl_vector_fixed_ref_const<T, num_cols> const& b) const
00595   {
00596     return vnl_matrix_fixed_mat_vec_mult<T,num_rows,num_cols>(*this,b);
00597   }
00598 #endif
00599 
00600   //: Set values of this matrix to those of M, starting at [top,left]
00601   vnl_matrix_fixed_ref const & update (vnl_matrix<T> const&, unsigned top=0, unsigned left=0) const;
00602 
00603   //: Set the elements of the i'th column to v[i]  (No bounds checking)
00604   vnl_matrix_fixed_ref const& set_column(unsigned i, T const * v) const;
00605 
00606   //: Set the elements of the i'th column to value, then return *this.
00607   vnl_matrix_fixed_ref const& set_column(unsigned i, T value ) const;
00608 
00609   //: Set j-th column to v, then return *this.
00610   vnl_matrix_fixed_ref const& set_column(unsigned j, vnl_vector<T> const& v) const;
00611 
00612   //: Set j-th column to v, then return *this.
00613   vnl_matrix_fixed_ref const& set_column(unsigned j, vnl_vector_fixed<T, num_rows> const& v) const;
00614 
00615   //: Set columns to those in M, starting at starting_column, then return *this.
00616   vnl_matrix_fixed_ref const& set_columns(unsigned starting_column, vnl_matrix<T> const& M) const;
00617 
00618   //: Set the elements of the i'th row to v[i]  (No bounds checking)
00619   vnl_matrix_fixed_ref const& set_row   (unsigned i, T const * v) const;
00620 
00621   //: Set the elements of the i'th row to value, then return *this.
00622   vnl_matrix_fixed_ref const& set_row   (unsigned i, T value ) const;
00623 
00624   //: Set the i-th row to v, then return *this.
00625   vnl_matrix_fixed_ref const& set_row   (unsigned i, vnl_vector<T> const& v) const;
00626 
00627   //: Set the i-th row to v, then return *this.
00628   vnl_matrix_fixed_ref const& set_row   (unsigned i, vnl_vector_fixed<T, num_cols> const& v) const;
00629 
00630   // ==== mutators ====
00631 
00632   //: Sets this matrix to an identity matrix, then returns "*this".
00633   //  Returning "*this" allows e.g. passing an identity matrix as argument to
00634   //  a function f, without having to name the constructed matrix:
00635   //  \code
00636   //     f(vnl_matrix_fixed_ref<double,5,5>().set_identity());
00637   //  \endcode
00638   //  Returning "*this" also allows "chaining" two or more operations:
00639   //  e.g., to set a 3x3 matrix to [3 0 0][0 2 0][0 0 1], one could say
00640   //  \code
00641   //     M.set_identity().scale_row(0,3).scale_column(1,2);
00642   //  \endcode
00643   //  If the matrix is not square, anyhow set main diagonal to 1, the rest to 0.
00644   vnl_matrix_fixed_ref const& set_identity() const;
00645 
00646   //: Reverses the order of rows, and returns "*this".
00647   //  Returning "*this" allows "chaining" two or more operations:
00648   //  e.g., to flip both up-down and left-right, one could just say
00649   //  \code
00650   //     M.flipud().fliplr();
00651   //  \endcode
00652   vnl_matrix_fixed_ref const& flipud() const;
00653 
00654   //: Reverses the order of columns, and returns "*this".
00655   //  Returning "*this" allows "chaining" two or more operations:
00656   //  e.g., to flip both up-down and left-right, one could just say
00657   //  \code
00658   //     M.flipud().fliplr();
00659   //  \endcode
00660   vnl_matrix_fixed_ref const& fliplr() const;
00661 
00662   //: Normalizes each row so it is a unit vector, and returns "*this".
00663   //  Zero rows are not modified
00664   //  Returning "*this" allows "chaining" two or more operations:
00665   //  e.g., to set a matrix to a row-normalized all-elements-equal matrix, say
00666   //  \code
00667   //     M.fill(1).normalize_rows();
00668   //  \endcode
00669   //  Returning "*this" also allows passing such a matrix as argument
00670   //  to a function f, without having to name the constructed matrix:
00671   //  \code
00672   //     f(vnl_matrix_fixed_ref<double,5,5>(1.0).normalize_rows());
00673   //  \endcode
00674   vnl_matrix_fixed_ref const& normalize_rows() const;
00675 
00676   //: Normalizes each column so it is a unit vector, and returns "*this".
00677   //  Zero columns are not modified
00678   //  Returning "*this" allows "chaining" two or more operations:
00679   //  e.g., to set a matrix to a column-normalized all-elements-equal matrix, say
00680   //  \code
00681   //     M.fill(1).normalize_columns();
00682   //  \endcode
00683   //  Returning "*this" also allows passing such a matrix as argument
00684   //  to a function f, without having to name the constructed matrix:
00685   //  \code
00686   //     f(vnl_matrix_fixed_ref<double,5,5>(1.0).normalize_columns());
00687   //  \endcode
00688   vnl_matrix_fixed_ref const& normalize_columns() const;
00689 
00690   //: Scales elements in given row by a factor T, and returns "*this".
00691   //  Returning "*this" allows "chaining" two or more operations:
00692   //  e.g., to set a 3x3 matrix to [3 0 0][0 2 0][0 0 1], one could say
00693   //  \code
00694   //     M.set_identity().scale_row(0,3).scale_column(1,2);
00695   //  \endcode
00696   vnl_matrix_fixed_ref const& scale_row   (unsigned row, T value) const;
00697 
00698   //: Scales elements in given column by a factor T, and returns "*this".
00699   //  Returning "*this" allows "chaining" two or more operations:
00700   //  e.g., to set a 3x3 matrix to [3 0 0][0 2 0][0 0 1], one could say
00701   //  \code
00702   //     M.set_identity().scale_row(0,3).scale_column(1,2);
00703   //  \endcode
00704   vnl_matrix_fixed_ref const& scale_column(unsigned col, T value) const;
00705 
00706   ////----------------------- Input/Output ----------------------------
00707 
00708   // : Read a vnl_matrix from an ascii vcl_istream, automatically determining file size if the input matrix has zero size.
00709   bool read_ascii(vcl_istream& s) const;
00710 
00711   //----------------------------------------------------------------------
00712   // Conversion to vnl_matrix_ref.
00713 
00714   // The const version of as_ref should return a const vnl_matrix_ref
00715   // so that the vnl_matrix_ref::non_const() cannot be used on
00716   // it. This prevents a const vnl_matrix_fixed from being cast into a
00717   // non-const vnl_matrix reference, giving a slight increase in type safety.
00718 
00719   //: Explicit conversion to a vnl_matrix_ref.
00720   // This is a cheap conversion for those functions that have an interface
00721   // for vnl_matrix_ref but not for vnl_matrix_fixed_ref. There is also a
00722   // conversion operator that should work most of the time.
00723   // \sa vnl_matrix_ref::non_const
00724   vnl_matrix_ref<T> as_ref() { return vnl_matrix_ref<T>( num_rows, num_cols, data_block() ); }
00725 
00726   //: Explicit conversion to a vnl_matrix_ref.
00727   // This is a cheap conversion for those functions that have an interface
00728   // for vnl_matrix_ref but not for vnl_matrix_fixed_ref. There is also a
00729   // conversion operator that should work most of the time.
00730   // \sa vnl_matrix_ref::non_const
00731   const vnl_matrix_ref<T> as_ref() const { return vnl_matrix_ref<T>( num_rows, num_cols, const_cast<T*>(data_block()) ); }
00732 
00733   //: Cheap conversion to vnl_matrix_ref
00734   // Sometimes, such as with templated functions, the compiler cannot
00735   // use this user-defined conversion. For those cases, use the
00736   // explicit as_ref() method instead.
00737   operator const vnl_matrix_ref<T>() const { return vnl_matrix_ref<T>( num_rows, num_cols, const_cast<T*>(data_block()) ); }
00738 
00739   //: Convert to a vnl_matrix.
00740   const vnl_matrix<T> as_matrix() const { return vnl_matrix<T>(const_cast<T*>(data_block()),num_rows,num_cols); }
00741 
00742   //----------------------------------------------------------------------
00743 
00744   typedef T element_type;
00745 
00746   //: Iterators
00747   typedef T       *iterator;
00748   //: Iterator pointing to start of data
00749   iterator       begin() const { return data_block(); }
00750   //: Iterator pointing to element beyond end of data
00751   iterator       end() const { return begin() + this->size(); }
00752   //--------------------------------------------------------------------------------
00753 
00754   //: Return true if *this == rhs
00755   bool operator_eq (vnl_matrix_fixed_ref_const<T,num_rows,num_cols> const & rhs) const
00756   {
00757     return vnl_matrix_fixed_ref<T,num_rows,num_cols>::equal( this->data_block(), rhs.data_block() );
00758   }
00759 
00760   //: Equality operator
00761   bool operator==(vnl_matrix_fixed_ref_const<T,num_rows,num_cols> const &that) const { return  this->operator_eq(that); }
00762 
00763   //: Inequality operator
00764   bool operator!=(vnl_matrix_fixed_ref_const<T,num_rows,num_cols> const &that) const { return !this->operator_eq(that); }
00765 
00766 //--------------------------------------------------------------------------------
00767 };
00768 
00769 #undef VNL_MATRIX_FIXED_VCL60_WORKAROUND
00770 
00771   // Helper routines for arithmetic. These routines know the size from
00772   // the template parameters. The vector-vector operations are
00773   // element-wise.
00774 
00775 // Make the operators below inline because (1) they are small and
00776 // (2) we then have less explicit instantiation trouble.
00777 
00778 // --- Matrix-scalar -------------------------------------------------------------
00779 
00780 template<class T, unsigned m, unsigned n>
00781 inline
00782 vnl_matrix_fixed<T,m,n> operator+( const vnl_matrix_fixed_ref_const<T,m,n>& mat1, const vnl_matrix_fixed_ref_const<T,m,n>& mat2 )
00783 {
00784   vnl_matrix_fixed<T,m,n> r;
00785   vnl_matrix_fixed<T,m,n>::add( mat1.data_block(), mat2.data_block(), r.data_block() );
00786   return r;
00787 }
00788 
00789 template<class T, unsigned m, unsigned n>
00790 inline
00791 vnl_matrix_fixed<T,m,n> operator+( const vnl_matrix_fixed_ref_const<T,m,n>& mat, T s )
00792 {
00793   vnl_matrix_fixed<T,m,n> r;
00794   vnl_matrix_fixed<T,m,n>::add( mat.data_block(), s, r.data_block() );
00795   return r;
00796 }
00797 
00798 template<class T, unsigned m, unsigned n>
00799 inline
00800 vnl_matrix_fixed<T,m,n> operator+( T s, const vnl_matrix_fixed_ref_const<T,m,n>& mat )
00801 {
00802   vnl_matrix_fixed<T,m,n> r;
00803   vnl_matrix_fixed<T,m,n>::add( mat.data_block(), s, r.data_block() );
00804   return r;
00805 }
00806 
00807 template<class T, unsigned m, unsigned n>
00808 inline
00809 vnl_matrix_fixed<T,m,n> operator-( const vnl_matrix_fixed_ref_const<T,m,n>& mat1, const vnl_matrix_fixed_ref_const<T,m,n>& mat2 )
00810 {
00811   vnl_matrix_fixed<T,m,n> r;
00812   vnl_matrix_fixed<T,m,n>::sub( mat1.data_block(), mat2.data_block(), r.data_block() );
00813   return r;
00814 }
00815 
00816 template<class T, unsigned m, unsigned n>
00817 inline
00818 vnl_matrix_fixed<T,m,n> operator-( const vnl_matrix_fixed_ref_const<T,m,n>& mat, T s )
00819 {
00820   vnl_matrix_fixed<T,m,n> r;
00821   vnl_matrix_fixed<T,m,n>::sub( mat.data_block(), s, r.data_block() );
00822   return r;
00823 }
00824 
00825 template<class T, unsigned m, unsigned n>
00826 inline
00827 vnl_matrix_fixed<T,m,n> operator-( T s, const vnl_matrix_fixed_ref_const<T,m,n>& mat )
00828 {
00829   vnl_matrix_fixed<T,m,n> r;
00830   vnl_matrix_fixed<T,m,n>::sub( s, mat.data_block(), r.data_block() );
00831   return r;
00832 }
00833 
00834 template<class T, unsigned m, unsigned n>
00835 inline
00836 vnl_matrix_fixed<T,m,n> operator*( const vnl_matrix_fixed_ref_const<T,m,n>& mat, T s )
00837 {
00838   vnl_matrix_fixed<T,m,n> r;
00839   vnl_matrix_fixed<T,m,n>::mul( mat.data_block(), s, r.data_block() );
00840   return r;
00841 }
00842 
00843 template<class T, unsigned m, unsigned n>
00844 inline
00845 vnl_matrix_fixed<T,m,n> operator*( T s, const vnl_matrix_fixed_ref_const<T,m,n>& mat )
00846 {
00847   vnl_matrix_fixed<T,m,n> r;
00848   vnl_matrix_fixed<T,m,n>::mul( mat.data_block(), s, r.data_block() );
00849   return r;
00850 }
00851 
00852 template<class T, unsigned m, unsigned n>
00853 inline
00854 vnl_matrix_fixed<T,m,n> operator/( const vnl_matrix_fixed_ref_const<T,m,n>& mat, T s )
00855 {
00856   vnl_matrix_fixed<T,m,n> r;
00857   vnl_matrix_fixed<T,m,n>::div( mat.data_block(), s, r.data_block() );
00858   return r;
00859 }
00860 
00861 
00862 template<class T, unsigned m, unsigned n>
00863 inline
00864 vnl_matrix_fixed<T,m,n> element_product( const vnl_matrix_fixed_ref_const<T,m,n>& mat1,
00865                                          const vnl_matrix_fixed_ref_const<T,m,n>& mat2 )
00866 {
00867   vnl_matrix_fixed<T,m,n> r;
00868   vnl_matrix_fixed<T,m,n>::mul( mat1.data_block(), mat2.data_block(), r.data_block() );
00869   return r;
00870 }
00871 
00872 
00873 template<class T, unsigned m, unsigned n>
00874 inline
00875 vnl_matrix_fixed<T,m,n> element_quotient( const vnl_matrix_fixed_ref_const<T,m,n>& mat1,
00876                                           const vnl_matrix_fixed_ref_const<T,m,n>& mat2)
00877 {
00878   vnl_matrix_fixed<T,m,n> r;
00879   vnl_matrix_fixed<T,m,n>::div( mat1.data_block(), mat2.data_block(), r.data_block() );
00880   return r;
00881 }
00882 
00883 
00884 // The following two functions are helper functions that keep the
00885 // matrix-matrix and matrix-vector multiplication code in one place,
00886 // so that bug fixes, etc, can be localized.
00887 template <class T, unsigned M, unsigned N>
00888 inline
00889 vnl_vector_fixed<T, M>
00890 vnl_matrix_fixed_mat_vec_mult(const vnl_matrix_fixed_ref_const<T, M, N>& a,
00891                               const vnl_vector_fixed_ref_const<T, N>& b)
00892 {
00893   vnl_vector_fixed<T, M> out;
00894   for (unsigned i = 0; i < M; ++i)
00895   {
00896     T accum = a(i,0) * b(0);
00897     for (unsigned k = 1; k < N; ++k)
00898       accum += a(i,k) * b(k);
00899     out(i) = accum;
00900   }
00901   return out;
00902 }
00903 
00904 // see comment above
00905 template <class T, unsigned M, unsigned N, unsigned O>
00906 inline
00907 vnl_matrix_fixed<T, M, O>
00908 vnl_matrix_fixed_mat_mat_mult(const vnl_matrix_fixed_ref_const<T, M, N>& a,
00909                               const vnl_matrix_fixed_ref_const<T, N, O>& b)
00910 {
00911   vnl_matrix_fixed<T, M, O> out;
00912   for (unsigned i = 0; i < M; ++i)
00913     for (unsigned j = 0; j < O; ++j)
00914     {
00915       T accum = a(i,0) * b(0,j);
00916       for (unsigned k = 1; k < N; ++k)
00917         accum += a(i,k) * b(k,j);
00918       out(i,j) = accum;
00919     }
00920   return out;
00921 }
00922 
00923 #ifndef VCL_VC_6
00924 // The version for correct compilers
00925 
00926 //: Multiply  conformant vnl_matrix_fixed (M x N) and vector_fixed (N)
00927 // \relatesalso vnl_vector_fixed
00928 // \relatesalso vnl_matrix_fixed
00929 template <class T, unsigned M, unsigned N>
00930 inline
00931 vnl_vector_fixed<T, M> operator*(const vnl_matrix_fixed_ref_const<T, M, N>& a, const vnl_vector_fixed_ref_const<T, N>& b)
00932 {
00933   return vnl_matrix_fixed_mat_vec_mult(a,b);
00934 }
00935 
00936 //: Multiply two conformant vnl_matrix_fixed (M x N) times (N x O)
00937 // \relatesalso vnl_matrix_fixed
00938 template <class T, unsigned M, unsigned N, unsigned O>
00939 inline
00940 vnl_matrix_fixed<T, M, O> operator*(const vnl_matrix_fixed_ref_const<T, M, N>& a, const vnl_matrix_fixed_ref_const<T, N, O>& b)
00941 {
00942   return vnl_matrix_fixed_mat_mat_mult(a,b);
00943 }
00944 #endif // ! VCL_VC_6
00945 
00946 
00947 // These overloads for the common case of mixing a fixed with a
00948 // non-fixed. Because the operator* are templated, the fixed will not
00949 // be automatically converted to a non-fixed-ref. These do it for you.
00950 
00951 template<class T, unsigned m, unsigned n>
00952 inline vnl_matrix<T> operator+( const vnl_matrix_fixed_ref_const<T,m,n>& a, const vnl_matrix<T>& b )
00953 {
00954   return a.as_ref() + b;
00955 }
00956 
00957 template<class T, unsigned m, unsigned n>
00958 inline vnl_matrix<T> operator+( const vnl_matrix<T>& a, const vnl_matrix_fixed_ref_const<T,m,n>& b )
00959 {
00960   return a + b.as_ref();
00961 }
00962 
00963 template<class T, unsigned m, unsigned n>
00964 inline vnl_matrix<T> operator-( const vnl_matrix_fixed_ref_const<T,m,n>& a, const vnl_matrix<T>& b )
00965 {
00966   return a.as_ref() - b;
00967 }
00968 
00969 template<class T, unsigned m, unsigned n>
00970 inline vnl_matrix<T> operator-( const vnl_matrix<T>& a, const vnl_matrix_fixed_ref_const<T,m,n>& b )
00971 {
00972   return a - b.as_ref();
00973 }
00974 
00975 template<class T, unsigned m, unsigned n>
00976 inline vnl_matrix<T> operator*( const vnl_matrix_fixed_ref_const<T,m,n>& a, const vnl_matrix<T>& b )
00977 {
00978   return a.as_ref() * b;
00979 }
00980 
00981 template<class T, unsigned m, unsigned n>
00982 inline vnl_matrix<T> operator*( const vnl_matrix<T>& a, const vnl_matrix_fixed_ref_const<T,m,n>& b )
00983 {
00984   return a * b.as_ref();
00985 }
00986 
00987 template<class T, unsigned m, unsigned n>
00988 inline vnl_vector<T> operator*( const vnl_matrix_fixed_ref_const<T,m,n>& a, const vnl_vector<T>& b )
00989 {
00990   return a.as_ref() * b;
00991 }
00992 
00993 template<class T, unsigned n>
00994 inline vnl_vector<T> operator*( const vnl_matrix<T>& a, const vnl_vector_fixed_ref_const<T,n>& b )
00995 {
00996   return a * b.as_ref();
00997 }
00998 
00999 
01000 // --- I/O operations ------------------------------------------------------------
01001 
01002 template<class T, unsigned m, unsigned n>
01003 inline
01004 vcl_ostream& operator<< (vcl_ostream& os, vnl_matrix_fixed_ref_const<T,m,n> const& mat)
01005 {
01006   mat.print(os);
01007   return os;
01008 }
01009 
01010 template<class T, unsigned m, unsigned n>
01011 inline
01012 vcl_istream& operator>> (vcl_istream& is, vnl_matrix_fixed_ref<T,m,n> const& mat)
01013 {
01014   mat.read_ascii(is);
01015   return is;
01016 }
01017 
01018 
01019 #endif // vnl_matrix_fixed_ref_h_