core/vnl/vnl_matrix_fixed.txx
Go to the documentation of this file.
00001 // This is core/vnl/vnl_matrix_fixed.txx
00002 #ifndef vnl_matrix_fixed_txx_
00003 #define vnl_matrix_fixed_txx_
00004 //:
00005 // \file
00006 #include "vnl_matrix_fixed.h"
00007 
00008 #include <vcl_cmath.h>
00009 #include <vcl_iostream.h>
00010 #include <vcl_cstdlib.h> // for abort
00011 #include <vcl_cassert.h>
00012 
00013 #include <vnl/vnl_error.h>
00014 #include <vnl/vnl_math.h>
00015 #include <vnl/vnl_vector_fixed.h>
00016 
00017 template<class T, unsigned nrows, unsigned ncols>
00018 void
00019 vnl_matrix_fixed<T,nrows,ncols>::add( const T* a, const T* b, T* r )
00020 {
00021   unsigned int count = nrows*ncols;
00022   while ( count-- )
00023     *(r++) = *(a++) + *(b++);
00024 }
00025 
00026 
00027 template<class T, unsigned nrows, unsigned ncols>
00028 void
00029 vnl_matrix_fixed<T,nrows,ncols>::add( const T* a, T b, T* r )
00030 {
00031   unsigned int count = nrows*ncols;
00032   while ( count-- )
00033     *(r++) = *(a++) + b;
00034 }
00035 
00036 template<class T, unsigned nrows, unsigned ncols>
00037 void
00038 vnl_matrix_fixed<T,nrows,ncols>::sub( const T* a, const T* b, T* r )
00039 {
00040   unsigned int count = nrows*ncols;
00041   while ( count-- )
00042     *(r++) = *(a++) - *(b++);
00043 }
00044 
00045 template<class T, unsigned nrows, unsigned ncols>
00046 void
00047 vnl_matrix_fixed<T,nrows,ncols>::sub( const T* a, T b, T* r )
00048 {
00049   unsigned int count = nrows*ncols;
00050   while ( count-- )
00051     *(r++) = *(a++) - b;
00052 }
00053 
00054 template<class T, unsigned nrows, unsigned ncols>
00055 void
00056 vnl_matrix_fixed<T,nrows,ncols>::sub( T a, const T* b, T* r )
00057 {
00058   unsigned int count = nrows*ncols;
00059   while ( count-- )
00060     *(r++) = a - *(b++);
00061 }
00062 
00063 template<class T, unsigned nrows, unsigned ncols>
00064 void
00065 vnl_matrix_fixed<T,nrows,ncols>::mul( const T* a, const T* b, T* r )
00066 {
00067   unsigned int count = nrows*ncols;
00068   while ( count-- )
00069     *(r++) = *(a++) * *(b++);
00070 }
00071 
00072 template<class T, unsigned nrows, unsigned ncols>
00073 void
00074 vnl_matrix_fixed<T,nrows,ncols>::mul( const T* a, T b, T* r )
00075 {
00076   unsigned int count = nrows*ncols;
00077   while ( count-- )
00078     *(r++) = *(a++) * b;
00079 }
00080 
00081 template<class T, unsigned nrows, unsigned ncols>
00082 void
00083 vnl_matrix_fixed<T,nrows,ncols>::div( const T* a, const T* b, T* r )
00084 {
00085   unsigned int count = nrows*ncols;
00086   while ( count-- )
00087     *(r++) = *(a++) / *(b++);
00088 }
00089 
00090 template<class T, unsigned nrows, unsigned ncols>
00091 void
00092 vnl_matrix_fixed<T,nrows,ncols>::div( const T* a, T b, T* r )
00093 {
00094   unsigned int count = nrows*ncols;
00095   while ( count-- )
00096     *(r++) = *(a++) / b;
00097 }
00098 
00099 template<class T, unsigned nrows, unsigned ncols>
00100 bool
00101 vnl_matrix_fixed<T,nrows,ncols>::equal( const T* a, const T* b )
00102 {
00103   unsigned int count = nrows*ncols;
00104   while ( count-- )
00105     if ( *(a++) != *(b++) )  return false;
00106   return true;
00107 }
00108 
00109 //------------------------------------------------------------
00110 
00111 
00112 template<class T, unsigned nrows, unsigned ncols>
00113 vnl_matrix_fixed<T,nrows,ncols>&
00114 vnl_matrix_fixed<T,nrows,ncols>::fill (T value)
00115 {
00116   for (unsigned int i = 0; i < nrows; ++i)
00117     for (unsigned int j = 0; j < ncols; ++j)
00118       this->data_[i][j] = value;
00119   return *this;
00120 }
00121 
00122 
00123 template<class T, unsigned nrows, unsigned ncols>
00124 vnl_matrix_fixed<T,nrows,ncols>&
00125 vnl_matrix_fixed<T,nrows,ncols>::fill_diagonal (T value)
00126 {
00127   for (unsigned int i = 0; i < nrows && i < ncols; ++i)
00128     this->data_[i][i] = value;
00129   return *this;
00130 }
00131 
00132 
00133 template<class T, unsigned nrows, unsigned ncols>
00134 vnl_matrix_fixed<T,nrows,ncols>&
00135 vnl_matrix_fixed<T,nrows,ncols>::set_diagonal(vnl_vector<T> const& diag)
00136 {
00137   assert(diag.size() >= nrows || diag.size() >= ncols);
00138   // The length of the diagonal of a non-square matrix is the minimum of
00139   // the matrix's width & height; that explains the "||" in the assert,
00140   // and the "&&" in the upper bound for the "for".
00141   for (unsigned int i = 0; i < nrows && i < ncols; ++i)
00142     this->data_[i][i] = diag[i];
00143   return *this;
00144 }
00145 
00146 
00147 template<class T, unsigned nrows, unsigned ncols>
00148 void
00149 vnl_matrix_fixed<T,nrows,ncols>::print(vcl_ostream& os) const
00150 {
00151   for (unsigned int i = 0; i < nrows; ++i)
00152   {
00153     os << this->data_[i][0];
00154     for (unsigned int j = 1; j < ncols; ++j)
00155       os << ' ' << this->data_[i][j];
00156     os << '\n';
00157   }
00158 }
00159 
00160 
00161 template <class T, unsigned nrows, unsigned ncols>
00162 vnl_matrix_fixed<T,nrows,ncols>
00163 vnl_matrix_fixed<T,nrows,ncols>::apply(T (*f)(T const&)) const
00164 {
00165   vnl_matrix_fixed<T,nrows,ncols> ret;
00166   vnl_c_vector<T>::apply(this->data_[0], rows()*cols(), f, ret.data_block());
00167   return ret;
00168 }
00169 
00170 template <class T, unsigned nrows, unsigned ncols>
00171 vnl_matrix_fixed<T,nrows,ncols>
00172 vnl_matrix_fixed<T,nrows,ncols>::apply(T (*f)(T)) const
00173 {
00174   vnl_matrix_fixed<T,nrows,ncols> ret;
00175   vnl_c_vector<T>::apply(this->data_[0], rows()*cols(), f, ret.data_block());
00176   return ret;
00177 }
00178 
00179 ////--------------------------- Additions------------------------------------
00180 
00181 
00182 template<class T, unsigned nrows, unsigned ncols>
00183 vnl_matrix_fixed<T,ncols,nrows>
00184 vnl_matrix_fixed<T,nrows,ncols>::transpose() const
00185 {
00186   vnl_matrix_fixed<T,ncols,nrows> result;
00187   for (unsigned int i = 0; i < cols(); ++i)
00188     for (unsigned int j = 0; j < rows(); ++j)
00189       result(i,j) = this->data_[j][i];
00190   return result;
00191 }
00192 
00193 template<class T, unsigned nrows, unsigned ncols>
00194 vnl_matrix_fixed<T,ncols,nrows>
00195 vnl_matrix_fixed<T,nrows,ncols>::conjugate_transpose() const
00196 {
00197   vnl_matrix_fixed<T,ncols,nrows> result(transpose());
00198   vnl_c_vector<T>::conjugate(result.begin(),  // src
00199                              result.begin(),  // dst
00200                              result.size());  // size of block
00201   return result;
00202 }
00203 
00204 template<class T, unsigned nrows, unsigned ncols>
00205 vnl_matrix_fixed<T,nrows,ncols>&
00206 vnl_matrix_fixed<T,nrows,ncols>::update (vnl_matrix<T> const& m,
00207                                          unsigned top, unsigned left)
00208 {
00209   const unsigned int bottom = top + m.rows();
00210   const unsigned int right = left + m.cols();
00211 #ifndef NDEBUG
00212   if (nrows < bottom || ncols < right)
00213     vnl_error_matrix_dimension ("update",
00214                                 bottom, right, m.rows(), m.cols());
00215 #endif
00216   for (unsigned int i = top; i < bottom; ++i)
00217     for (unsigned int j = left; j < right; ++j)
00218       this->data_[i][j] = m(i-top,j-left);
00219   return *this;
00220 }
00221 
00222 
00223 template<class T, unsigned nrows, unsigned ncols>
00224 vnl_matrix<T>
00225 vnl_matrix_fixed<T,nrows,ncols>::extract (unsigned rowz, unsigned colz,
00226                                           unsigned top, unsigned left) const
00227 {
00228   vnl_matrix<T> result(rowz, colz);
00229   this->extract( result, top, left );
00230   return result;
00231 }
00232 
00233 
00234 template<class T, unsigned nrows, unsigned ncols>
00235 void
00236 vnl_matrix_fixed<T,nrows,ncols>::extract (vnl_matrix<T>& sub_matrix,
00237                                           unsigned top, unsigned left) const
00238 {
00239   unsigned int rowz = sub_matrix.rows();
00240   unsigned int colz = sub_matrix.cols();
00241 #ifndef NDEBUG
00242   unsigned int bottom = top + rowz;
00243   unsigned int right = left + colz;
00244   if ((nrows < bottom) || (ncols < right))
00245     vnl_error_matrix_dimension ("extract",
00246                                 nrows, ncols, bottom, right);
00247 #endif
00248   for (unsigned int i = 0; i < rowz; ++i)      // actual copy of all elements
00249     for (unsigned int j = 0; j < colz; ++j)    // in submatrix
00250       sub_matrix(i,j) = this->data_[top+i][left+j];
00251 }
00252 
00253 
00254 template<class T, unsigned nrows, unsigned ncols>
00255 vnl_matrix_fixed<T,nrows,ncols>&
00256 vnl_matrix_fixed<T,nrows,ncols>::copy_in(T const *p)
00257 {
00258   T* dp = this->data_block();
00259   unsigned int i = nrows * ncols;
00260   while (i--)
00261     *dp++ = *p++;
00262   return *this;
00263 }
00264 
00265 template<class T, unsigned nrows, unsigned ncols>
00266 void vnl_matrix_fixed<T,nrows,ncols>::copy_out(T *p) const
00267 {
00268   T const* dp = this->data_block();
00269   unsigned int i = nrows*ncols;
00270   while (i--)
00271     *p++ = *dp++;
00272 }
00273 
00274 template<class T, unsigned nrows, unsigned ncols>
00275 vnl_matrix_fixed<T,nrows,ncols>&
00276 vnl_matrix_fixed<T,nrows,ncols>::set_identity()
00277 {
00278   // Two simple loops are generally better than having a branch inside
00279   // the loop. Probably worth the O(n) extra writes.
00280   for (unsigned int i = 0; i < nrows; ++i)
00281     for (unsigned int j = 0; j < ncols; ++j)
00282       this->data_[i][j] = T(0);
00283   for (unsigned int i = 0; i < nrows && i < ncols; ++i)
00284     this->data_[i][i] = T(1);
00285   return *this;
00286 }
00287 
00288 //: Make each row of the matrix have unit norm.
00289 // All-zero rows are ignored.
00290 template<class T, unsigned nrows, unsigned ncols>
00291 vnl_matrix_fixed<T,nrows,ncols>&
00292 vnl_matrix_fixed<T,nrows,ncols>::normalize_rows()
00293 {
00294   for (unsigned int i = 0; i < nrows; ++i)
00295   {
00296     abs_t norm(0); // double will not do for all types.
00297     for (unsigned int j = 0; j < ncols; ++j)
00298       norm += vnl_math_squared_magnitude( this->data_[i][j] );
00299 
00300     if (norm != 0)
00301     {
00302       typedef typename vnl_numeric_traits<abs_t>::real_t real_t;
00303       real_t scale = real_t(1)/vcl_sqrt((real_t)norm);
00304       for (unsigned int j = 0; j < ncols; ++j)
00305       {
00306         // FIXME need correct rounding here
00307         // There is e.g. no *standard* operator*=(complex<float>, double), hence the T() cast.
00308         this->data_[i][j] *= T(scale);
00309       }
00310     }
00311   }
00312   return *this;
00313 }
00314 
00315 template<class T, unsigned nrows, unsigned ncols>
00316 vnl_matrix_fixed<T,nrows,ncols>&
00317 vnl_matrix_fixed<T,nrows,ncols>::normalize_columns()
00318 {
00319   for (unsigned int j = 0; j < ncols; ++j) {  // For each column in the Matrix
00320     abs_t norm(0); // double will not do for all types.
00321     for (unsigned int i = 0; i < nrows; ++i)
00322       norm += vnl_math_squared_magnitude( this->data_[i][j] );
00323 
00324     if (norm != 0)
00325     {
00326       typedef typename vnl_numeric_traits<abs_t>::real_t real_t;
00327       real_t scale = real_t(1)/vcl_sqrt((real_t)norm);
00328       for (unsigned int i = 0; i < nrows; ++i)
00329       {
00330         // FIXME need correct rounding here
00331         // There is e.g. no *standard* operator*=(complex<float>, double), hence the T() cast.
00332         this->data_[i][j] *= T(scale);
00333       }
00334     }
00335   }
00336   return *this;
00337 }
00338 
00339 template<class T, unsigned nrows, unsigned ncols>
00340 vnl_matrix_fixed<T,nrows,ncols>&
00341 vnl_matrix_fixed<T,nrows,ncols>::scale_row(unsigned row_index, T value)
00342 {
00343 #ifndef NDEBUG
00344   if (row_index >= nrows)
00345     vnl_error_matrix_row_index("scale_row", row_index);
00346 #endif
00347   for (unsigned int j = 0; j < ncols; ++j)
00348     this->data_[row_index][j] *= value;
00349   return *this;
00350 }
00351 
00352 template<class T, unsigned nrows, unsigned ncols>
00353 vnl_matrix_fixed<T,nrows,ncols>&
00354 vnl_matrix_fixed<T,nrows,ncols>::scale_column(unsigned column_index, T value)
00355 {
00356 #ifndef NDEBUG
00357   if (column_index >= ncols)
00358     vnl_error_matrix_col_index("scale_column", column_index);
00359 #endif
00360   for (unsigned int j = 0; j < nrows; ++j)
00361     this->data_[j][column_index] *= value;
00362   return *this;
00363 }
00364 
00365 //: Returns a copy of n rows, starting from "row"
00366 template<class T, unsigned nrows, unsigned ncols>
00367 vnl_matrix<T>
00368 vnl_matrix_fixed<T,nrows,ncols>::get_n_rows (unsigned row, unsigned n) const
00369 {
00370 #ifndef NDEBUG
00371   if (row + n > nrows)
00372     vnl_error_matrix_row_index ("get_n_rows", row);
00373 #endif
00374 
00375   // Extract data rowwise.
00376   return vnl_matrix<T>(data_[row], n, ncols);
00377 }
00378 
00379 template<class T, unsigned nrows, unsigned ncols>
00380 vnl_matrix<T>
00381 vnl_matrix_fixed<T,nrows,ncols>::get_n_columns (unsigned column, unsigned n) const
00382 {
00383 #ifndef NDEBUG
00384   if (column + n > ncols)
00385     vnl_error_matrix_col_index ("get_n_columns", column);
00386 #endif
00387 
00388   vnl_matrix<T> result(nrows, n);
00389   for (unsigned int c = 0; c < n; ++c)
00390     for (unsigned int r = 0; r < nrows; ++r)
00391       result(r, c) = this->data_[r][column + c];
00392   return result;
00393 }
00394 
00395 //: Create a vector out of row[row_index].
00396 template<class T, unsigned nrows, unsigned ncols>
00397 vnl_vector_fixed<T,ncols> vnl_matrix_fixed<T,nrows,ncols>::get_row(unsigned row_index) const
00398 {
00399 #ifdef ERROR_CHECKING
00400   if (row_index >= nrows)
00401     vnl_error_matrix_row_index ("get_row", row_index);
00402 #endif
00403 
00404   vnl_vector_fixed<T,ncols> v;
00405   for (unsigned int j = 0; j < ncols; ++j)    // For each element in row
00406     v[j] = this->data_[row_index][j];
00407   return v;
00408 }
00409 
00410 //: Create a vector out of column[column_index].
00411 template<class T, unsigned nrows, unsigned ncols>
00412 vnl_vector_fixed<T,nrows> vnl_matrix_fixed<T,nrows,ncols>::get_column(unsigned column_index) const
00413 {
00414 #ifdef ERROR_CHECKING
00415   if (column_index >= ncols)
00416     vnl_error_matrix_col_index ("get_column", column_index);
00417 #endif
00418 
00419   vnl_vector_fixed<T,nrows> v;
00420   for (unsigned int j = 0; j < nrows; ++j)
00421     v[j] = this->data_[j][column_index];
00422   return v;
00423 }
00424 
00425 //: Return a vector with the content of the (main) diagonal
00426 template<class T, unsigned nrows, unsigned ncols>
00427 vnl_vector<T> vnl_matrix_fixed<T,nrows,ncols>::get_diagonal() const
00428 {
00429   vnl_vector<T> v(nrows < ncols ? nrows : ncols);
00430   for (unsigned int j = 0; j < nrows && j < ncols; ++j)
00431     v[j] = this->data_[j][j];
00432   return v;
00433 }
00434 
00435 //--------------------------------------------------------------------------------
00436 
00437 template<class T, unsigned nrows, unsigned ncols>
00438 vnl_matrix_fixed<T,nrows,ncols>&
00439 vnl_matrix_fixed<T,nrows,ncols>::set_row(unsigned row_index, T const *v)
00440 {
00441   for (unsigned int j = 0; j < ncols; ++j)
00442     this->data_[row_index][j] = v[j];
00443   return *this;
00444 }
00445 
00446 template<class T, unsigned nrows, unsigned ncols>
00447 vnl_matrix_fixed<T,nrows,ncols>&
00448 vnl_matrix_fixed<T,nrows,ncols>::set_row(unsigned row_index, vnl_vector<T> const &v)
00449 {
00450   if (v.size() >= ncols)
00451     set_row(row_index,v.data_block());
00452   else
00453     for (unsigned int j = 0; j < v.size(); ++j)
00454       this->data_[row_index][j] = v[j];
00455   return *this;
00456 }
00457 
00458 template<class T, unsigned nrows, unsigned ncols>
00459 vnl_matrix_fixed<T,nrows,ncols>&
00460 vnl_matrix_fixed<T,nrows,ncols>::set_row(unsigned row_index, vnl_vector_fixed<T,ncols> const &v)
00461 {
00462   set_row(row_index,v.data_block());
00463   return *this;
00464 }
00465 
00466 template<class T, unsigned nrows, unsigned ncols>
00467 vnl_matrix_fixed<T,nrows,ncols>&
00468 vnl_matrix_fixed<T,nrows,ncols>::set_row(unsigned row_index, T v)
00469 {
00470   for (unsigned int j = 0; j < ncols; ++j)
00471     this->data_[row_index][j] = v;
00472   return *this;
00473 }
00474 
00475 //--------------------------------------------------------------------------------
00476 
00477 template<class T, unsigned nrows, unsigned ncols>
00478 vnl_matrix_fixed<T,nrows,ncols>&
00479 vnl_matrix_fixed<T,nrows,ncols>::set_column(unsigned column_index, T const *v)
00480 {
00481   for (unsigned int i = 0; i < nrows; ++i)
00482     this->data_[i][column_index] = v[i];
00483   return *this;
00484 }
00485 
00486 template<class T, unsigned nrows, unsigned ncols>
00487 vnl_matrix_fixed<T,nrows,ncols>&
00488 vnl_matrix_fixed<T,nrows,ncols>::set_column(unsigned column_index, vnl_vector<T> const &v)
00489 {
00490   if (v.size() >= nrows)
00491     set_column(column_index,v.data_block());
00492   else
00493     for (unsigned int i = 0; i < v.size(); ++i)
00494       this->data_[i][column_index] = v[i];
00495   return *this;
00496 }
00497 
00498 template<class T, unsigned nrows, unsigned ncols>
00499 vnl_matrix_fixed<T,nrows,ncols>&
00500 vnl_matrix_fixed<T,nrows,ncols>::set_column(unsigned column_index, vnl_vector_fixed<T,nrows> const &v)
00501 {
00502   set_column(column_index,v.data_block());
00503   return *this;
00504 }
00505 
00506 template<class T, unsigned nrows, unsigned ncols>
00507 vnl_matrix_fixed<T,nrows,ncols>&
00508 vnl_matrix_fixed<T,nrows,ncols>::set_column(unsigned column_index, T v)
00509 {
00510   for (unsigned int j = 0; j < nrows; ++j)
00511     this->data_[j][column_index] = v;
00512   return *this;
00513 }
00514 
00515 
00516 template<class T, unsigned nrows, unsigned ncols>
00517 vnl_matrix_fixed<T,nrows,ncols>&
00518 vnl_matrix_fixed<T,nrows,ncols>::set_columns(unsigned starting_column, vnl_matrix<T> const& m)
00519 {
00520   for (unsigned int j = 0; j < m.cols() && starting_column+j < ncols; ++j) // don't go too far right; possibly only use part of m
00521     for (unsigned int i = 0; i < nrows && i < m.rows(); ++i) // smallest of the two heights; possibly only use part of m
00522       this->data_[i][starting_column + j] = m(i,j);
00523   return *this;
00524 }
00525 
00526 
00527 template <class T, unsigned nrows, unsigned ncols>
00528 bool
00529 vnl_matrix_fixed<T,nrows,ncols>::is_identity() const
00530 {
00531   T const zero(0);
00532   T const one(1);
00533   for (unsigned int i = 0; i < nrows; ++i)
00534     for (unsigned int j = 0; j < ncols; ++j)
00535     {
00536       T xm = this->data_[i][j];
00537       if ( !((i == j) ? (xm == one) : (xm == zero)) )
00538         return false;
00539     }
00540   return true;
00541 }
00542 
00543 //: Return true if maximum absolute deviation of M from identity is <= tol.
00544 template <class T, unsigned nrows, unsigned ncols>
00545 bool
00546 vnl_matrix_fixed<T,nrows,ncols>::is_identity(double tol) const
00547 {
00548   T one(1);
00549   for (unsigned int i = 0; i < nrows; ++i)
00550     for (unsigned int j = 0; j < ncols; ++j)
00551     {
00552       T xm = this->data_[i][j];
00553       abs_t absdev = (i == j) ? vnl_math_abs(xm - one) : vnl_math_abs(xm);
00554       if (absdev > tol)
00555         return false;
00556     }
00557   return true;
00558 }
00559 
00560 template <class T, unsigned nrows, unsigned ncols>
00561 bool
00562 vnl_matrix_fixed<T,nrows,ncols>::is_zero() const
00563 {
00564   T const zero(0);
00565   for (unsigned int i = 0; i < nrows; ++i)
00566     for (unsigned int j = 0; j < ncols; ++j)
00567       if ( !( this->data_[i][ j] == zero) )
00568         return false;
00569 
00570   return true;
00571 }
00572 
00573 template <class T, unsigned nrows, unsigned ncols>
00574 bool
00575 vnl_matrix_fixed<T,nrows,ncols>::is_zero(double tol) const
00576 {
00577   for (unsigned int i = 0; i < nrows; ++i)
00578     for (unsigned int j = 0; j < ncols; ++j)
00579       if (vnl_math_abs(this->data_[i][j]) > tol)
00580         return false;
00581 
00582   return true;
00583 }
00584 
00585 template <class T, unsigned nrows, unsigned ncols>
00586 bool
00587 vnl_matrix_fixed<T,nrows,ncols>::has_nans() const
00588 {
00589   for (unsigned int i = 0; i < nrows; ++i)
00590     for (unsigned int j = 0; j < ncols; ++j)
00591       if (vnl_math_isnan(this->data_[i][j]))
00592         return true;
00593 
00594   return false;
00595 }
00596 
00597 template <class T, unsigned nrows, unsigned ncols>
00598 bool
00599 vnl_matrix_fixed<T,nrows,ncols>::is_finite() const
00600 {
00601   for (unsigned int i = 0; i < nrows; ++i)
00602     for (unsigned int j = 0; j < ncols; ++j)
00603       if (!vnl_math_isfinite(this->data_[i][j]))
00604         return false;
00605 
00606   return true;
00607 }
00608 
00609 //: Abort if any element of M is inf or nan
00610 template <class T, unsigned nrows, unsigned ncols>
00611 void
00612 vnl_matrix_fixed<T,nrows,ncols>::assert_finite_internal() const
00613 {
00614   if (is_finite())
00615     return;
00616 
00617   vcl_cerr << "\n\n" __FILE__ ": " << __LINE__ << ": matrix has non-finite elements\n";
00618 
00619   if (rows() <= 20 && cols() <= 20)
00620     vcl_cerr << __FILE__ ": here it is:\n" << *this << '\n';
00621   else
00622   {
00623     vcl_cerr << __FILE__ ": it is quite big (" << rows() << 'x' << cols() << ")\n"
00624              << __FILE__ ": in the following picture '-' means finite and '*' means non-finite:\n";
00625 
00626     for (unsigned int i=0; i<rows(); ++i)
00627     {
00628       for (unsigned int j=0; j<cols(); ++j)
00629         vcl_cerr << char(vnl_math_isfinite(this->data_[i][ j]) ? '-' : '*');
00630       vcl_cerr << '\n';
00631     }
00632   }
00633   vcl_cerr << __FILE__ ": calling abort()\n";
00634   vcl_abort();
00635 }
00636 
00637 //: Abort unless M has the given size.
00638 template <class T, unsigned nrows, unsigned ncols>
00639 void
00640 vnl_matrix_fixed<T,nrows,ncols>::assert_size_internal(unsigned rs,unsigned cs) const
00641 {
00642   if (nrows!=rs || ncols!=cs)
00643   {
00644     vcl_cerr << __FILE__ ": size is " << nrows << 'x' << ncols
00645              << ". should be " << rs << 'x' << cs << vcl_endl;
00646     vcl_abort();
00647   }
00648 }
00649 
00650 template <class T, unsigned nrows, unsigned ncols>
00651 bool
00652 vnl_matrix_fixed<T,nrows,ncols>::read_ascii(vcl_istream& s)
00653 {
00654   if (!s.good())
00655   {
00656     vcl_cerr << __FILE__ ": vnl_matrix_fixed<T,nrows,ncols>::read_ascii: Called with bad stream\n";
00657     return false;
00658   }
00659 
00660   for (unsigned int i = 0; i < nrows; ++i)
00661     for (unsigned int j = 0; j < ncols; ++j)
00662       s >> this->data_[i][j];
00663 
00664   return s.good() || s.eof();
00665 }
00666 
00667 
00668 template <class T, unsigned nrows, unsigned ncols>
00669 vnl_matrix_fixed<T,nrows,ncols>&
00670 vnl_matrix_fixed<T,nrows,ncols>::flipud()
00671 {
00672   for (unsigned int r1 = 0; 2*r1+1 < nrows; ++r1)
00673   {
00674     unsigned int r2 = nrows - 1 - r1;
00675     for (unsigned int c = 0; c < ncols; ++c)
00676     {
00677       T tmp = this->data_[r1][c];
00678       this->data_[r1][c] = this->data_[r2][c];
00679       this->data_[r2][c] = tmp;
00680     }
00681   }
00682   return *this;
00683 }
00684 
00685 
00686 template <class T, unsigned nrows, unsigned ncols>
00687 vnl_matrix_fixed<T,nrows,ncols>&
00688 vnl_matrix_fixed<T,nrows,ncols>::fliplr()
00689 {
00690   for (unsigned int c1 = 0; 2*c1+1 < ncols; ++c1)
00691   {
00692     unsigned int c2 = ncols - 1 - c1;
00693     for (unsigned int r = 0; r < nrows; ++r)
00694     {
00695       T tmp = this->data_[r][c1];
00696       this->data_[r][c1] = this->data_[r][c2];
00697       this->data_[r][c2] = tmp;
00698     }
00699   }
00700   return *this;
00701 }
00702 
00703 template <class T, unsigned nrows, unsigned ncols>
00704 typename vnl_matrix_fixed<T,nrows,ncols>::abs_t
00705 vnl_matrix_fixed<T,nrows,ncols>::operator_one_norm() const
00706 {
00707   abs_t m(0);
00708   for (unsigned int j=0; j<ncols; ++j)
00709   {
00710     abs_t t(0);
00711     for (unsigned int i=0; i<nrows; ++i)
00712       t += vnl_math_abs( this->data_[i][j] );
00713     if (t > m)
00714       m = t;
00715   }
00716   return m;
00717 }
00718 
00719 template <class T, unsigned nrows, unsigned ncols>
00720 typename vnl_matrix_fixed<T,nrows,ncols>::abs_t
00721 vnl_matrix_fixed<T,nrows,ncols>::operator_inf_norm() const
00722 {
00723   abs_t m(0);
00724   for (unsigned int i=0; i<nrows; ++i)
00725   {
00726     abs_t t(0);
00727     for (unsigned int j=0; j<ncols; ++j)
00728       t += vnl_math_abs( this->data_[i][j] );
00729     if (t > m)
00730       m = t;
00731   }
00732   return m;
00733 }
00734 
00735 //: Transpose square matrix M in place.
00736 template <class T, unsigned nrows, unsigned ncols>
00737 vnl_matrix_fixed<T,nrows,ncols>&
00738 vnl_matrix_fixed<T,nrows,ncols>::inplace_transpose()
00739 {
00740   assert(nrows==ncols); // cannot inplace_transpose non-square fixed size matrix
00741   for (unsigned i = 0; i < nrows; ++i)
00742   for (unsigned j = i+1; j < ncols; ++j)
00743   {
00744     T t = this->data_[i][j];
00745     this->data_[i][j] = this->data_[j][i];
00746     this->data_[j][i] = t;
00747   }
00748   return *this;
00749 }
00750 
00751 // Workaround for argument deduction bug in VC6. See comment in .h
00752 // file.  Note that the body of the function is outside the #ifdefs to
00753 // maintain a single implementation of the function. The way to read
00754 // this code is to just jump to the #else part and scan down, ignoring
00755 // the #endif. Unless, of course, you are masochist and actually want
00756 // to read the workaround.
00757 //
00758 #if defined(VCL_VC_6) && !defined(__GCCXML__)
00759 
00760 template<class VecA, class VecB, class RM>
00761 RM
00762 outer_product_fixed_calc_helper<VecA,VecB,RM>::calc( VecA const& a, VecB const& b )
00763 {
00764   RM out; // RM should be a vnl_matrix_fixed of VecA::SIZE by VecB::SIZE
00765   for (unsigned int i = 0; i < VecA::SIZE; ++i)
00766     for (unsigned int j = 0; j < VecB::SIZE; ++j)
00767       out[i][j] = a[i] * b[j];
00768   return out;
00769 };
00770 
00771 #define VNL_OUTER_PRODUCT_FIXED_INSTANTIATE( T, M, N ) \
00772 template struct outer_product_fixed_calc_helper<vnl_vector_fixed<T,M >, \
00773                                                 vnl_vector_fixed<T,N >, \
00774                                                 vnl_matrix_fixed<T,M,N > >
00775 
00776 #else // no need for workaround; declare the function sanely.
00777 
00778 template <class T, unsigned m, unsigned n>
00779 vnl_matrix_fixed<T,m,n>
00780 outer_product(vnl_vector_fixed<T,m> const& a, vnl_vector_fixed<T,n> const& b)
00781 {
00782   vnl_matrix_fixed<T,m,n> out; // = a.column() * b.row()
00783   for (unsigned int i = 0; i < m; ++i)
00784     for (unsigned int j = 0; j < n; ++j)
00785       out[i][j] = a[i] * b[j];
00786   return out;
00787 }
00788 
00789 #define VNL_OUTER_PRODUCT_FIXED_INSTANTIATE( T, M, N ) \
00790 template vnl_matrix_fixed<T,M,N > outer_product(vnl_vector_fixed<T,M > const&,\
00791                                                 vnl_vector_fixed<T,N > const& )
00792 
00793 #endif // VC60 outer_product workaround
00794 
00795 
00796 #undef VNL_MATRIX_FIXED_INSTANTIATE
00797 #define VNL_MATRIX_FIXED_INSTANTIATE(T, M, N) \
00798 template class vnl_matrix_fixed<T,M,N >; \
00799 VNL_OUTER_PRODUCT_FIXED_INSTANTIATE( T, M, N )
00800 
00801 #endif // vnl_matrix_fixed_txx_