00001 // This is core/vnl/vnl_matrix.h 00002 #ifndef vnl_matrix_h_ 00003 #define vnl_matrix_h_ 00004 #ifdef VCL_NEEDS_PRAGMA_INTERFACE 00005 #pragma interface 00006 #endif 00007 //: 00008 // \file 00009 // \brief An ordinary mathematical matrix 00010 // \verbatim 00011 // Modifications 00012 // Apr 21, 1989 - MBN - Initial design and implementation 00013 // Jun 22, 1989 - MBN - Removed non-destructive methods 00014 // Aug 09, 1989 - LGO - Inherit from Generic 00015 // Aug 20, 1989 - MBN - Changed template usage to reflect new syntax 00016 // Sep 11, 1989 - MBN - Added conditional exception handling and base class 00017 // Oct 05, 1989 - LGO - Don't re-allocate data in operator= when same size 00018 // Oct 19, 1989 - LGO - Add extra parameter to varargs constructor 00019 // Oct 19, 1989 - MBN - Added optional argument to set_compare method 00020 // Dec 08, 1989 - LGO - Allocate column data in one chunk 00021 // Dec 08, 1989 - LGO - Clean-up get and put, add const everywhere. 00022 // Dec 19, 1989 - LGO - Remove the map and reduce methods 00023 // Feb 22, 1990 - MBN - Changed size arguments from int to unsigned int 00024 // Jun 30, 1990 - MJF - Added base class name to constructor initializer 00025 // Feb 21, 1992 - VDN - New lite version 00026 // May 05, 1992 - VDN - Use envelope to avoid unnecessary copying 00027 // Sep 30, 1992 - VDN - Matrix inversion with singular value decomposition 00028 // Aug 21, 1996 - AWF - set_identity, normalize_rows, scale_row. 00029 // Sep 30, 1996 - AWF - set_row/column methods. Const-correct data_block(). 00030 // 14 Feb 1997 - AWF - get_n_rows, get_n_columns. 00031 // 20 Mar 1997 - PVR - get_row, get_column. 00032 // 24-Oct-2010 - Peter Vanroose - mutators and filling methods now return *this 00033 // 18-Jan-2011 - Peter Vanroose - added methods set_diagonal() & get_diagonal() 00034 // \endverbatim 00035 00036 #include <vcl_iosfwd.h> 00037 #include <vnl/vnl_tag.h> 00038 #include <vnl/vnl_c_vector.h> 00039 #include <vnl/vnl_config.h> 00040 #ifndef NDEBUG 00041 # if VNL_CONFIG_CHECK_BOUNDS 00042 # include <vnl/vnl_error.h> 00043 # include <vcl_cassert.h> 00044 # endif 00045 #else 00046 # undef VNL_CONFIG_CHECK_BOUNDS 00047 # define VNL_CONFIG_CHECK_BOUNDS 0 00048 # undef ERROR_CHECKING 00049 #endif 00050 00051 export template <class T> class vnl_vector; 00052 export template <class T> class vnl_matrix; 00053 00054 //-------------------------------------------------------------------------------- 00055 00056 #ifndef DOXYGEN_SHOULD_SKIP_THIS 00057 #define v vnl_vector<T> 00058 #define m vnl_matrix<T> 00059 #endif // DOXYGEN_SHOULD_SKIP_THIS 00060 template <class T> m operator+(T const&, m const&); 00061 template <class T> m operator-(T const&, m const&); 00062 template <class T> m operator*(T const&, m const&); 00063 template <class T> m element_product(m const&, m const&); 00064 template <class T> m element_quotient(m const&, m const&); 00065 template <class T> T dot_product(m const&, m const&); 00066 template <class T> T inner_product(m const&, m const&); 00067 template <class T> T cos_angle(m const&, m const& ); 00068 template <class T> vcl_ostream& operator<<(vcl_ostream&, m const&); 00069 template <class T> vcl_istream& operator>>(vcl_istream&, m&); 00070 #undef v 00071 #undef m 00072 00073 //-------------------------------------------------------------------------------- 00074 00075 enum vnl_matrix_type 00076 { 00077 vnl_matrix_null, 00078 vnl_matrix_identity 00079 }; 00080 00081 //: An ordinary mathematical matrix 00082 // The vnl_matrix<T> class implements two-dimensional arithmetic 00083 // matrices for a user-specified numeric data type. Using the 00084 // parameterized types facility of C++, it is possible, for 00085 // example, for the user to create a matrix of rational numbers 00086 // by parameterizing the vnl_matrix class over the Rational class. 00087 // The only requirement for the type is that it supports the 00088 // basic arithmetic operators. 00089 // 00090 // Note: Unlike the other sequence classes, the 00091 // vnl_matrix<T> class is fixed-size. It will not grow once the 00092 // size has been specified to the constructor or changed by the 00093 // assignment or multiplication operators. The vnl_matrix<T> 00094 // class is row-based with addresses of rows being cached, and 00095 // elements accessed as m[row][col]. 00096 // 00097 // Note: The matrix can, however, be resized using the set_size(nr,nc) function. 00098 // 00099 // Note: Indexing of the matrix is zero-based, so the top-left element is M(0,0). 00100 // 00101 // Note: Inversion of matrix M, and other operations such as solving systems of linear 00102 // equations are handled by the matrix decomposition classes in vnl/algo, such 00103 // as matrix_inverse, svd, qr etc. 00104 // 00105 // Note: Use a vnl_vector<T> with these matrices. 00106 00107 template<class T> 00108 class vnl_matrix 00109 { 00110 public: 00111 //: Default constructor creates an empty matrix of size 0,0. 00112 vnl_matrix() : 00113 num_rows(0), 00114 num_cols(0), 00115 data(0) 00116 { 00117 } 00118 00119 //: Construct a matrix of size r rows by c columns 00120 // Contents are unspecified. 00121 // Complexity $O(1)$ 00122 vnl_matrix(unsigned r, unsigned c); // r rows, c cols. 00123 00124 //: Construct a matrix of size r rows by c columns, and all elements equal to v0 00125 // Complexity $O(r.c)$ 00126 vnl_matrix(unsigned r, unsigned c, T const& v0); // r rows, c cols, value v0. 00127 00128 //: Construct a matrix of size r rows by c columns, with a special type 00129 // Contents are specified by t 00130 // Complexity $O(r.c)$ 00131 vnl_matrix(unsigned r, unsigned c, vnl_matrix_type t); // r rows, c cols, special type 00132 00133 //: Construct a matrix of size r rows by c columns, initialised by an automatic array 00134 // The first n elements, are initialised row-wise, to values. 00135 // Complexity $O(n)$ 00136 vnl_matrix(unsigned r, unsigned c, unsigned n, T const values[]); // use automatic arrays. 00137 00138 //: Construct a matrix of size r rows by c columns, initialised by a memory block 00139 // The values are initialise row wise from the data. 00140 // Complexity $O(r.c)$ 00141 vnl_matrix(T const* data_block, unsigned r, unsigned c); // fill row-wise. 00142 00143 //: Copy construct a matrix 00144 // Complexity $O(r.c)$ 00145 vnl_matrix(vnl_matrix<T> const&); // from another matrix. 00146 00147 #ifndef VXL_DOXYGEN_SHOULD_SKIP_THIS 00148 // <internal> 00149 // These constructors are here so that operator* etc can take 00150 // advantage of the C++ return value optimization. 00151 vnl_matrix(vnl_matrix<T> const &, vnl_matrix<T> const &, vnl_tag_add); // M + M 00152 vnl_matrix(vnl_matrix<T> const &, vnl_matrix<T> const &, vnl_tag_sub); // M - M 00153 vnl_matrix(vnl_matrix<T> const &, T, vnl_tag_mul); // M * s 00154 vnl_matrix(vnl_matrix<T> const &, T, vnl_tag_div); // M / s 00155 vnl_matrix(vnl_matrix<T> const &, T, vnl_tag_add); // M + s 00156 vnl_matrix(vnl_matrix<T> const &, T, vnl_tag_sub); // M - s 00157 vnl_matrix(vnl_matrix<T> const &, vnl_matrix<T> const &, vnl_tag_mul); // M * M 00158 vnl_matrix(vnl_matrix<T> &that, vnl_tag_grab) 00159 : num_rows(that.num_rows), num_cols(that.num_cols), data(that.data) 00160 { that.num_cols=that.num_rows=0; that.data=0; } // "*this" now uses "that"'s data. 00161 // </internal> 00162 #endif 00163 00164 //: Matrix destructor 00165 ~vnl_matrix(); 00166 00167 // Basic 2D-Array functionality------------------------------------------- 00168 00169 //: Return number of rows 00170 unsigned rows() const { return num_rows; } 00171 00172 //: Return number of columns 00173 // A synonym for cols() 00174 unsigned columns() const { return num_cols; } 00175 00176 //: Return number of columns 00177 // A synonym for columns() 00178 unsigned cols() const { return num_cols; } 00179 00180 //: Return number of elements 00181 // This equals rows() * cols() 00182 unsigned size() const { return rows()*cols(); } 00183 00184 //: set element with boundary checks if error checking is on. 00185 void put(unsigned r, unsigned c, T const&); 00186 00187 //: get element with boundary checks if error checking is on. 00188 T get(unsigned r, unsigned c) const; 00189 00190 //: return pointer to given row 00191 // No boundary checking here. 00192 T * operator[](unsigned r) { return data[r]; } 00193 00194 //: return pointer to given row 00195 // No boundary checking here. 00196 T const * operator[](unsigned r) const { return data[r]; } 00197 00198 //: Access an element for reading or writing 00199 // There are assert style boundary checks - #define NDEBUG to turn them off. 00200 T & operator()(unsigned r, unsigned c) 00201 { 00202 #if VNL_CONFIG_CHECK_BOUNDS 00203 assert(r<rows()); // Check the row index is valid 00204 assert(c<cols()); // Check the column index is valid 00205 #endif 00206 return this->data[r][c]; 00207 } 00208 00209 //: Access an element for reading 00210 // There are assert style boundary checks - #define NDEBUG to turn them off. 00211 T const & operator()(unsigned r, unsigned c) const 00212 { 00213 #if VNL_CONFIG_CHECK_BOUNDS 00214 assert(r<rows()); // Check the row index is valid 00215 assert(c<cols()); // Check the column index is valid 00216 #endif 00217 return this->data[r][c]; 00218 } 00219 00220 00221 // ----------------------- Filling and copying ----------------------- 00222 00223 //: Sets all elements of matrix to specified value, and returns "*this". 00224 // Complexity $O(r.c)$ 00225 // Returning "*this" allows "chaining" two or more operations: 00226 // e.g., to set a matrix to a column-normalized all-elements-equal matrix, say 00227 // \code 00228 // M.fill(1).normalize_columns(); 00229 // \endcode 00230 // Returning "*this" also allows passing such a matrix as argument 00231 // to a function f, without having to name the constructed matrix: 00232 // \code 00233 // f(vnl_matrix<double>(5,5,1.0).normalize_columns()); 00234 // \endcode 00235 vnl_matrix& fill(T const&); 00236 00237 //: Sets all diagonal elements of matrix to specified value; returns "*this". 00238 // Complexity $O(\min(r,c))$ 00239 // Returning "*this" allows "chaining" two or more operations: 00240 // e.g., to set a 3x3 matrix to [5 0 0][0 10 0][0 0 15], just say 00241 // \code 00242 // M.fill_diagonal(5).scale_row(1,2).scale_column(2,3); 00243 // \endcode 00244 // Returning "*this" also allows passing a diagonal-filled matrix as argument 00245 // to a function f, without having to name the constructed matrix: 00246 // \code 00247 // f(vnl_matrix<double>(3,3).fill_diagonal(5)); 00248 // \endcode 00249 vnl_matrix& fill_diagonal(T const&); 00250 00251 //: Sets the diagonal elements of this matrix to the specified list of values. 00252 // Returning "*this" allows "chaining" two or more operations: see the 00253 // reasoning (and the examples) in the documentation for method 00254 // fill_diagonal(). 00255 vnl_matrix& set_diagonal(vnl_vector<T> const&); 00256 00257 //: Fills (laminates) this matrix with the given data, then returns it. 00258 // We assume that the argument points to a contiguous rows*cols array, stored rowwise. 00259 // No bounds checking on the array. 00260 // Returning "*this" allows "chaining" two or more operations: 00261 // e.g., to fill a square matrix column-wise, fill it rowwise then transpose: 00262 // \code 00263 // M.copy_in(array).inplace_transpose(); 00264 // \endcode 00265 // Returning "*this" also allows passing a filled-in matrix as argument 00266 // to a function f, without having to name the constructed matrix: 00267 // \code 00268 // f(vnl_matrix<double>(3,3).copy_in(array)); 00269 // \endcode 00270 vnl_matrix& copy_in(T const *); 00271 00272 //: Fills (laminates) this matrix with the given data, then returns it. 00273 // A synonym for copy_in() 00274 vnl_matrix& set(T const *d) { return copy_in(d); } 00275 00276 //: Fills the given array with this matrix. 00277 // We assume that the argument points to a contiguous rows*cols array, stored rowwise. 00278 // No bounds checking on the array. 00279 void copy_out(T *) const; 00280 00281 //: Set all elements to value v 00282 // Complexity $O(r.c)$ 00283 vnl_matrix<T>& operator=(T const&v) { fill(v); return *this; } 00284 00285 //: Copies all elements of rhs matrix into lhs matrix. 00286 // Complexity $O(\min(r,c))$ 00287 vnl_matrix<T>& operator=(vnl_matrix<T> const&); 00288 00289 // ----------------------- Arithmetic -------------------------------- 00290 // note that these functions should not pass scalar as a const&. 00291 // Look what would happen to A /= A(0,0). 00292 00293 //: Add rhs to each element of lhs matrix in situ 00294 vnl_matrix<T>& operator+=(T value); 00295 00296 //: Subtract rhs from each element of lhs matrix in situ 00297 vnl_matrix<T>& operator-=(T value); 00298 00299 //: Scalar multiplication in situ of lhs matrix by rhs 00300 vnl_matrix<T>& operator*=(T value); 00301 00302 //: Scalar division of lhs matrix in situ by rhs 00303 vnl_matrix<T>& operator/=(T value); 00304 00305 //: Add rhs to lhs matrix in situ 00306 vnl_matrix<T>& operator+=(vnl_matrix<T> const&); 00307 //: Subtract rhs from lhs matrix in situ 00308 vnl_matrix<T>& operator-=(vnl_matrix<T> const&); 00309 //: Multiply lhs matrix in situ by rhs 00310 vnl_matrix<T>& operator*=(vnl_matrix<T> const&rhs) { return *this = (*this) * rhs; } 00311 00312 //: Negate all elements of matrix 00313 vnl_matrix<T> operator-() const; 00314 00315 00316 //: Add rhs to each element of lhs matrix and return result in new matrix 00317 vnl_matrix<T> operator+(T const& v) const { return vnl_matrix<T>(*this, v, vnl_tag_add()); } 00318 00319 //: Subtract rhs from each element of lhs matrix and return result in new matrix 00320 vnl_matrix<T> operator-(T const& v) const { return vnl_matrix<T>(*this, v, vnl_tag_sub()); } 00321 00322 //: Scalar multiplication of lhs matrix by rhs and return result in new matrix 00323 vnl_matrix<T> operator*(T const& v) const { return vnl_matrix<T>(*this, v, vnl_tag_mul()); } 00324 00325 //: Scalar division of lhs matrix by rhs and return result in new matrix 00326 vnl_matrix<T> operator/(T const& v) const { return vnl_matrix<T>(*this, v, vnl_tag_div()); } 00327 00328 //: Matrix add rhs to lhs matrix and return result in new matrix 00329 vnl_matrix<T> operator+(vnl_matrix<T> const& rhs) const { return vnl_matrix<T>(*this, rhs, vnl_tag_add()); } 00330 //: Matrix subtract rhs from lhs and return result in new matrix 00331 vnl_matrix<T> operator-(vnl_matrix<T> const& rhs) const { return vnl_matrix<T>(*this, rhs, vnl_tag_sub()); } 00332 //: Matrix multiply lhs by rhs matrix and return result in new matrix 00333 vnl_matrix<T> operator*(vnl_matrix<T> const& rhs) const { return vnl_matrix<T>(*this, rhs, vnl_tag_mul()); } 00334 00335 ////--------------------------- Additions ---------------------------- 00336 00337 //: Make a new matrix by applying function to each element. 00338 vnl_matrix<T> apply(T (*f)(T)) const; 00339 00340 //: Make a new matrix by applying function to each element. 00341 vnl_matrix<T> apply(T (*f)(T const&)) const; 00342 00343 //: Return transpose 00344 vnl_matrix<T> transpose() const; 00345 00346 //: Return conjugate transpose 00347 vnl_matrix<T> conjugate_transpose() const; 00348 00349 //: Set values of this matrix to those of M, starting at [top,left] 00350 vnl_matrix<T>& update(vnl_matrix<T> const&, unsigned top=0, unsigned left=0); 00351 00352 //: Set the elements of the i'th column to v[i] (No bounds checking) 00353 vnl_matrix& set_column(unsigned i, T const * v); 00354 00355 //: Set the elements of the i'th column to value, then return *this. 00356 vnl_matrix& set_column(unsigned i, T value ); 00357 00358 //: Set j-th column to v, then return *this. 00359 vnl_matrix& set_column(unsigned j, vnl_vector<T> const& v); 00360 00361 //: Set columns to those in M, starting at starting_column, then return *this. 00362 vnl_matrix& set_columns(unsigned starting_column, vnl_matrix<T> const& M); 00363 00364 //: Set the elements of the i'th row to v[i] (No bounds checking) 00365 vnl_matrix& set_row(unsigned i, T const * v); 00366 00367 //: Set the elements of the i'th row to value, then return *this. 00368 vnl_matrix& set_row(unsigned i, T value ); 00369 00370 //: Set the i-th row 00371 vnl_matrix& set_row(unsigned i, vnl_vector<T> const&); 00372 00373 //: Extract a sub-matrix of size r x c, starting at (top,left) 00374 // Thus it contains elements [top,top+r-1][left,left+c-1] 00375 vnl_matrix<T> extract(unsigned r, unsigned c, 00376 unsigned top=0, unsigned left=0) const; 00377 00378 //: Extract a sub-matrix starting at (top,left) 00379 // 00380 // The output is stored in \a sub_matrix, and it should have the 00381 // required size on entry. Thus the result will contain elements 00382 // [top,top+sub_matrix.rows()-1][left,left+sub_matrix.cols()-1] 00383 void extract ( vnl_matrix<T>& sub_matrix, 00384 unsigned top=0, unsigned left=0) const; 00385 00386 00387 //: Get a vector equal to the given row 00388 vnl_vector<T> get_row(unsigned r) const; 00389 00390 //: Get a vector equal to the given column 00391 vnl_vector<T> get_column(unsigned c) const; 00392 00393 //: Get n rows beginning at rowstart 00394 vnl_matrix<T> get_n_rows(unsigned rowstart, unsigned n) const; 00395 00396 //: Get n columns beginning at colstart 00397 vnl_matrix<T> get_n_columns(unsigned colstart, unsigned n) const; 00398 00399 //: Return a vector with the content of the (main) diagonal 00400 vnl_vector<T> get_diagonal() const; 00401 00402 // ==== mutators ==== 00403 00404 //: Sets this matrix to an identity matrix, then returns "*this". 00405 // Returning "*this" allows e.g. passing an identity matrix as argument to 00406 // a function f, without having to name the constructed matrix: 00407 // \code 00408 // f(vnl_matrix<double>(5,5).set_identity()); 00409 // \endcode 00410 // Returning "*this" also allows "chaining" two or more operations: 00411 // e.g., to set a 3x3 matrix to [3 0 0][0 2 0][0 0 1], one could say 00412 // \code 00413 // M.set_identity().scale_row(0,3).scale_column(1,2); 00414 // \endcode 00415 // If the matrix is not square, anyhow set main diagonal to 1, the rest to 0. 00416 vnl_matrix& set_identity(); 00417 00418 //: Transposes this matrix efficiently, and returns it. 00419 // Returning "*this" allows "chaining" two or more operations: 00420 // e.g., to fill a square matrix column-wise, fill it rowwise then transpose: 00421 // \code 00422 // M.copy_in(array).inplace_transpose(); 00423 // \endcode 00424 vnl_matrix& inplace_transpose(); 00425 00426 //: Reverses the order of rows, and returns "*this". 00427 // Returning "*this" allows "chaining" two or more operations: 00428 // e.g., to flip both up-down and left-right, one could just say 00429 // \code 00430 // M.flipud().fliplr(); 00431 // \endcode 00432 vnl_matrix& flipud(); 00433 00434 //: Reverses the order of columns, and returns "*this". 00435 // Returning "*this" allows "chaining" two or more operations: 00436 // e.g., to flip both up-down and left-right, one could just say 00437 // \code 00438 // M.flipud().fliplr(); 00439 // \endcode 00440 vnl_matrix& fliplr(); 00441 00442 //: Normalizes each row so it is a unit vector, and returns "*this". 00443 // Zero rows are not modified 00444 // Returning "*this" allows "chaining" two or more operations: 00445 // e.g., to set a matrix to a row-normalized all-elements-equal matrix, say 00446 // \code 00447 // M.fill(1).normalize_rows(); 00448 // \endcode 00449 // Returning "*this" also allows passing such a matrix as argument 00450 // to a function f, without having to name the constructed matrix: 00451 // \code 00452 // f(vnl_matrix<double>(5,5,1.0).normalize_rows()); 00453 // \endcode 00454 vnl_matrix& normalize_rows(); 00455 00456 //: Normalizes each column so it is a unit vector, and returns "*this". 00457 // Zero columns are not modified 00458 // Returning "*this" allows "chaining" two or more operations: 00459 // e.g., to set a matrix to a column-normalized all-elements-equal matrix, say 00460 // \code 00461 // M.fill(1).normalize_columns(); 00462 // \endcode 00463 // Returning "*this" also allows passing such a matrix as argument 00464 // to a function f, without having to name the constructed matrix: 00465 // \code 00466 // f(vnl_matrix<double>(5,5,1.0).normalize_columns()); 00467 // \endcode 00468 vnl_matrix& normalize_columns(); 00469 00470 //: Scales elements in given row by a factor T, and returns "*this". 00471 // Returning "*this" allows "chaining" two or more operations: 00472 // e.g., to set a 3x3 matrix to [3 0 0][0 2 0][0 0 1], one could say 00473 // \code 00474 // M.set_identity().scale_row(0,3).scale_column(1,2); 00475 // \endcode 00476 vnl_matrix& scale_row(unsigned row, T value); 00477 00478 //: Scales elements in given column by a factor T, and returns "*this". 00479 // Returning "*this" allows "chaining" two or more operations: 00480 // e.g., to set a 3x3 matrix to [3 0 0][0 2 0][0 0 1], one could say 00481 // \code 00482 // M.set_identity().scale_row(0,3).scale_column(1,2); 00483 // \endcode 00484 vnl_matrix& scale_column(unsigned col, T value); 00485 00486 //: Swap this matrix with that matrix 00487 void swap(vnl_matrix<T> & that); 00488 00489 //: Type def for norms. 00490 typedef typename vnl_c_vector<T>::abs_t abs_t; 00491 00492 //: Return sum of absolute values of elements 00493 abs_t array_one_norm() const { return vnl_c_vector<T>::one_norm(begin(), size()); } 00494 00495 //: Return square root of sum of squared absolute element values 00496 abs_t array_two_norm() const { return vnl_c_vector<T>::two_norm(begin(), size()); } 00497 00498 //: Return largest absolute element value 00499 abs_t array_inf_norm() const { return vnl_c_vector<T>::inf_norm(begin(), size()); } 00500 00501 //: Return sum of absolute values of elements 00502 abs_t absolute_value_sum() const { return array_one_norm(); } 00503 00504 //: Return largest absolute value 00505 abs_t absolute_value_max() const { return array_inf_norm(); } 00506 00507 // $ || M ||_1 := \max_j \sum_i | M_{ij} | $ 00508 abs_t operator_one_norm() const; 00509 00510 // $ || M ||_\inf := \max_i \sum_j | M_{ij} | $ 00511 abs_t operator_inf_norm() const; 00512 00513 //: Return Frobenius norm of matrix (sqrt of sum of squares of its elements) 00514 abs_t frobenius_norm() const { return vnl_c_vector<T>::two_norm(begin(), size()); } 00515 00516 //: Return Frobenius norm of matrix (sqrt of sum of squares of its elements) 00517 abs_t fro_norm() const { return frobenius_norm(); } 00518 00519 //: Return RMS of all elements 00520 abs_t rms() const { return vnl_c_vector<T>::rms_norm(begin(), size()); } 00521 00522 //: Return minimum value of elements 00523 T min_value() const { return vnl_c_vector<T>::min_value(begin(), size()); } 00524 00525 //: Return maximum value of elements 00526 T max_value() const { return vnl_c_vector<T>::max_value(begin(), size()); } 00527 00528 //: Return location of minimum value of elements 00529 unsigned arg_min() const { return vnl_c_vector<T>::arg_min(begin(), size()); } 00530 00531 //: Return location of maximum value of elements 00532 unsigned arg_max() const { return vnl_c_vector<T>::arg_max(begin(), size()); } 00533 00534 //: Return mean of all matrix elements 00535 T mean() const { return vnl_c_vector<T>::mean(begin(), size()); } 00536 00537 // predicates 00538 00539 //: Return true iff the size is zero. 00540 bool empty() const { return !data || !num_rows || !num_cols; } 00541 00542 //: Return true if all elements equal to identity. 00543 bool is_identity() const; 00544 00545 //: Return true if all elements equal to identity, within given tolerance 00546 bool is_identity(double tol) const; 00547 00548 //: Return true if all elements equal to zero. 00549 bool is_zero() const; 00550 00551 //: Return true if all elements equal to zero, within given tolerance 00552 bool is_zero(double tol) const; 00553 00554 //: Return true if all elements of both matrices are equal, within given tolerance 00555 bool is_equal(vnl_matrix<T> const& rhs, double tol) const; 00556 00557 //: Return true if finite 00558 bool is_finite() const; 00559 00560 //: Return true if matrix contains NaNs 00561 bool has_nans() const; 00562 00563 //: abort if size is not as expected 00564 // This function does or tests nothing if NDEBUG is defined 00565 void assert_size(unsigned r, unsigned c) const 00566 { 00567 #ifndef NDEBUG 00568 assert_size_internal(r, c); 00569 #endif 00570 } 00571 //: abort if matrix contains any INFs or NANs. 00572 // This function does or tests nothing if NDEBUG is defined 00573 void assert_finite() const 00574 { 00575 #ifndef NDEBUG 00576 assert_finite_internal(); 00577 #endif 00578 } 00579 00580 ////----------------------- Input/Output ---------------------------- 00581 00582 //: Read a vnl_matrix from an ascii vcl_istream, automatically determining file size if the input matrix has zero size. 00583 static vnl_matrix<T> read(vcl_istream& s); 00584 00585 // : Read a vnl_matrix from an ascii vcl_istream, automatically determining file size if the input matrix has zero size. 00586 bool read_ascii(vcl_istream& s); 00587 00588 //-------------------------------------------------------------------------------- 00589 00590 //: Access the contiguous block storing the elements in the matrix row-wise. O(1). 00591 // 1d array, row-major order. 00592 T const* data_block() const { return data[0]; } 00593 00594 //: Access the contiguous block storing the elements in the matrix row-wise. O(1). 00595 // 1d array, row-major order. 00596 T * data_block() { return data[0]; } 00597 00598 //: Access the 2D array, so that elements can be accessed with array[row][col] directly. 00599 // 2d array, [row][column]. 00600 T const* const* data_array() const { return data; } 00601 00602 //: Access the 2D array, so that elements can be accessed with array[row][col] directly. 00603 // 2d array, [row][column]. 00604 T * * data_array() { return data; } 00605 00606 typedef T element_type; 00607 00608 //: Iterators 00609 typedef T *iterator; 00610 //: Iterator pointing to start of data 00611 iterator begin() { return data?data[0]:0; } 00612 //: Iterator pointing to element beyond end of data 00613 iterator end() { return data?data[0]+num_rows*num_cols:0; } 00614 00615 //: Const iterators 00616 typedef T const *const_iterator; 00617 //: Iterator pointing to start of data 00618 const_iterator begin() const { return data?data[0]:0; } 00619 //: Iterator pointing to element beyond end of data 00620 const_iterator end() const { return data?data[0]+num_rows*num_cols:0; } 00621 00622 //: Return a reference to this. 00623 // Useful in code which would prefer not to know if its argument 00624 // is a matrix, matrix_ref or a matrix_fixed. Note that it doesn't 00625 // return a matrix_ref, so it's only useful in templates or macros. 00626 vnl_matrix<T> const& as_ref() const { return *this; } 00627 00628 //: Return a reference to this. 00629 vnl_matrix<T>& as_ref() { return *this; } 00630 00631 //-------------------------------------------------------------------------------- 00632 00633 //: Return true if *this == rhs 00634 bool operator_eq(vnl_matrix<T> const & rhs) const; 00635 00636 //: Equality operator 00637 bool operator==(vnl_matrix<T> const &that) const { return this->operator_eq(that); } 00638 00639 //: Inequality operator 00640 bool operator!=(vnl_matrix<T> const &that) const { return !this->operator_eq(that); } 00641 00642 //: Print matrix to os in some hopefully sensible format 00643 void print(vcl_ostream& os) const; 00644 00645 //: Make the matrix as if it had been default-constructed. 00646 void clear(); 00647 00648 //: Resize to r rows by c columns. Old data lost. 00649 // Returns true if size changed. 00650 bool set_size(unsigned r, unsigned c); 00651 00652 //-------------------------------------------------------------------------------- 00653 00654 protected: 00655 unsigned num_rows; // Number of rows 00656 unsigned num_cols; // Number of columns 00657 T** data; // Pointer to the vnl_matrix 00658 00659 #if VCL_HAS_SLICED_DESTRUCTOR_BUG 00660 // Since this bug exists, we need a flag that can be set during 00661 // construction to tell our destructor whether we own data. 00662 char vnl_matrix_own_data; 00663 #endif 00664 00665 void assert_size_internal(unsigned r, unsigned c) const; 00666 void assert_finite_internal() const; 00667 00668 //: Delete data 00669 void destroy(); 00670 00671 #if VCL_NEED_FRIEND_FOR_TEMPLATE_OVERLOAD 00672 # define v vnl_vector<T> 00673 # define m vnl_matrix<T> 00674 friend m operator+ VCL_NULL_TMPL_ARGS (T const&, m const&); 00675 friend m operator- VCL_NULL_TMPL_ARGS (T const&, m const&); 00676 friend m operator* VCL_NULL_TMPL_ARGS (T const&, m const&); 00677 friend m element_product VCL_NULL_TMPL_ARGS (m const&, m const&); 00678 friend m element_quotient VCL_NULL_TMPL_ARGS (m const&, m const&); 00679 friend T dot_product VCL_NULL_TMPL_ARGS (m const&, m const&); 00680 friend T inner_product VCL_NULL_TMPL_ARGS (m const&, m const&); 00681 friend T cos_angle VCL_NULL_TMPL_ARGS (m const&, m const&); 00682 friend vcl_ostream& operator<< VCL_NULL_TMPL_ARGS (vcl_ostream&, m const&); 00683 friend vcl_istream& operator>> VCL_NULL_TMPL_ARGS (vcl_istream&, m&); 00684 # undef v 00685 # undef m 00686 #endif 00687 00688 // inline function template instantiation hack for gcc 2.97 -- fsm 00689 static void inline_function_tickler(); 00690 }; 00691 00692 00693 // Definitions of inline functions. 00694 00695 00696 //: Returns the value of the element at specified row and column. O(1). 00697 // Checks for valid range of indices. 00698 00699 template<class T> 00700 inline T vnl_matrix<T>::get(unsigned row, unsigned column) const 00701 { 00702 #ifdef ERROR_CHECKING 00703 if (row >= this->num_rows) // If invalid size specified 00704 vnl_error_matrix_row_index("get", row); // Raise exception 00705 if (column >= this->num_cols) // If invalid size specified 00706 vnl_error_matrix_col_index("get", column); // Raise exception 00707 #endif 00708 return this->data[row][column]; 00709 } 00710 00711 //: Puts value into element at specified row and column. O(1). 00712 // Checks for valid range of indices. 00713 00714 template<class T> 00715 inline void vnl_matrix<T>::put(unsigned row, unsigned column, T const& value) 00716 { 00717 #ifdef ERROR_CHECKING 00718 if (row >= this->num_rows) // If invalid size specified 00719 vnl_error_matrix_row_index("put", row); // Raise exception 00720 if (column >= this->num_cols) // If invalid size specified 00721 vnl_error_matrix_col_index("put", column); // Raise exception 00722 #endif 00723 this->data[row][column] = value; // Assign data value 00724 } 00725 00726 00727 // non-member arithmetical operators. 00728 00729 //: 00730 // \relatesalso vnl_matrix 00731 template<class T> 00732 inline vnl_matrix<T> operator*(T const& value, vnl_matrix<T> const& m) 00733 { 00734 return vnl_matrix<T>(m, value, vnl_tag_mul()); 00735 } 00736 00737 //: 00738 // \relatesalso vnl_matrix 00739 template<class T> 00740 inline vnl_matrix<T> operator+(T const& value, vnl_matrix<T> const& m) 00741 { 00742 return vnl_matrix<T>(m, value, vnl_tag_add()); 00743 } 00744 00745 //: Swap two matrices 00746 // \relatesalso vnl_matrix 00747 template<class T> 00748 inline void swap(vnl_matrix<T> &A, vnl_matrix<T> &B) { A.swap(B); } 00749 00750 00751 #endif // vnl_matrix_h_