[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 2003-2008 by Gunnar Kedenburg and Ullrich Koethe */ 00004 /* */ 00005 /* This file is part of the VIGRA computer vision library. */ 00006 /* The VIGRA Website is */ 00007 /* http://hci.iwr.uni-heidelberg.de/vigra/ */ 00008 /* Please direct questions, bug reports, and contributions to */ 00009 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00010 /* vigra@informatik.uni-hamburg.de */ 00011 /* */ 00012 /* Permission is hereby granted, free of charge, to any person */ 00013 /* obtaining a copy of this software and associated documentation */ 00014 /* files (the "Software"), to deal in the Software without */ 00015 /* restriction, including without limitation the rights to use, */ 00016 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00017 /* sell copies of the Software, and to permit persons to whom the */ 00018 /* Software is furnished to do so, subject to the following */ 00019 /* conditions: */ 00020 /* */ 00021 /* The above copyright notice and this permission notice shall be */ 00022 /* included in all copies or substantial portions of the */ 00023 /* Software. */ 00024 /* */ 00025 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00026 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00027 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00028 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00029 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00030 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00031 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00032 /* OTHER DEALINGS IN THE SOFTWARE. */ 00033 /* */ 00034 /************************************************************************/ 00035 00036 00037 #ifndef VIGRA_MATRIX_HXX 00038 #define VIGRA_MATRIX_HXX 00039 00040 #include <cmath> 00041 #include <iosfwd> 00042 #include <iomanip> 00043 #include "multi_array.hxx" 00044 #include "mathutil.hxx" 00045 #include "numerictraits.hxx" 00046 00047 00048 namespace vigra 00049 { 00050 00051 /** \defgroup LinearAlgebraModule Linear Algebra 00052 00053 \brief Classes and functions for matrix algebra, linear equations systems, eigen systems, least squares etc. 00054 */ 00055 00056 /** \ingroup LinearAlgebraModule 00057 00058 Namespace <tt>vigra/linalg</tt> hold VIGRA's linear algebra functionality. But most of its contents 00059 is exported into namespace <tt>vigra</tt> via <tt>using</tt> directives. 00060 */ 00061 namespace linalg 00062 { 00063 00064 template <class T, class C> 00065 inline MultiArrayIndex 00066 rowCount(const MultiArrayView<2, T, C> &x); 00067 00068 template <class T, class C> 00069 inline MultiArrayIndex 00070 columnCount(const MultiArrayView<2, T, C> &x); 00071 00072 template <class T, class C> 00073 inline MultiArrayView <2, T, C> 00074 rowVector(MultiArrayView <2, T, C> const & m, MultiArrayIndex d); 00075 00076 template <class T, class C> 00077 inline MultiArrayView <2, T, C> 00078 columnVector(MultiArrayView<2, T, C> const & m, MultiArrayIndex d); 00079 00080 template <class T, class ALLOC> 00081 class TemporaryMatrix; 00082 00083 template <class T, class C1, class C2> 00084 void transpose(const MultiArrayView<2, T, C1> &v, MultiArrayView<2, T, C2> &r); 00085 00086 template <class T, class C> 00087 bool isSymmetric(const MultiArrayView<2, T, C> &v); 00088 00089 enum RawArrayMemoryLayout { RowMajor, ColumnMajor }; 00090 00091 /********************************************************/ 00092 /* */ 00093 /* Matrix */ 00094 /* */ 00095 /********************************************************/ 00096 00097 /** Matrix class. 00098 00099 \ingroup LinearAlgebraModule 00100 00101 This is the basic class for all linear algebra computations. Matrices are 00102 strored in a <i>column-major</i> format, i.e. the row index is varying fastest. 00103 This is the same format as in the lapack and gmm++ libraries, so it will 00104 be easy to interface these libraries. In fact, if you need optimized 00105 high performance code, you should use them. The VIGRA linear algebra 00106 functionality is provided for smaller problems and rapid prototyping 00107 (no one wants to spend half the day installing a new library just to 00108 discover that the new algorithm idea didn't work anyway). 00109 00110 <b>See also:</b> 00111 <ul> 00112 <li> \ref LinearAlgebraFunctions 00113 </ul> 00114 00115 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 00116 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00117 Namespaces: vigra and vigra::linalg 00118 */ 00119 template <class T, class ALLOC = std::allocator<T> > 00120 class Matrix 00121 : public MultiArray<2, T, ALLOC> 00122 { 00123 typedef MultiArray<2, T, ALLOC> BaseType; 00124 00125 public: 00126 typedef Matrix<T, ALLOC> matrix_type; 00127 typedef TemporaryMatrix<T, ALLOC> temp_type; 00128 typedef MultiArrayView<2, T, UnstridedArrayTag> view_type; 00129 typedef typename BaseType::value_type value_type; 00130 typedef typename BaseType::pointer pointer; 00131 typedef typename BaseType::const_pointer const_pointer; 00132 typedef typename BaseType::reference reference; 00133 typedef typename BaseType::const_reference const_reference; 00134 typedef typename BaseType::difference_type difference_type; 00135 typedef typename BaseType::difference_type_1 difference_type_1; 00136 typedef ALLOC allocator_type; 00137 00138 /** default constructor 00139 */ 00140 Matrix() 00141 {} 00142 00143 /** construct with given allocator 00144 */ 00145 explicit Matrix(ALLOC const & alloc) 00146 : BaseType(alloc) 00147 {} 00148 00149 /** construct with given shape and init all 00150 elements with zero. Note that the order of the axes is 00151 <tt>difference_type(rows, columns)</tt> which 00152 is the opposite of the usual VIGRA convention. 00153 */ 00154 explicit Matrix(const difference_type &shape, 00155 ALLOC const & alloc = allocator_type()) 00156 : BaseType(shape, alloc) 00157 {} 00158 00159 /** construct with given shape and init all 00160 elements with zero. Note that the order of the axes is 00161 <tt>(rows, columns)</tt> which 00162 is the opposite of the usual VIGRA convention. 00163 */ 00164 Matrix(difference_type_1 rows, difference_type_1 columns, 00165 ALLOC const & alloc = allocator_type()) 00166 : BaseType(difference_type(rows, columns), alloc) 00167 {} 00168 00169 /** construct with given shape and init all 00170 elements with the constant \a init. Note that the order of the axes is 00171 <tt>difference_type(rows, columns)</tt> which 00172 is the opposite of the usual VIGRA convention. 00173 */ 00174 Matrix(const difference_type &shape, const_reference init, 00175 allocator_type const & alloc = allocator_type()) 00176 : BaseType(shape, init, alloc) 00177 {} 00178 00179 /** construct with given shape and init all 00180 elements with the constant \a init. Note that the order of the axes is 00181 <tt>(rows, columns)</tt> which 00182 is the opposite of the usual VIGRA convention. 00183 */ 00184 Matrix(difference_type_1 rows, difference_type_1 columns, const_reference init, 00185 allocator_type const & alloc = allocator_type()) 00186 : BaseType(difference_type(rows, columns), init, alloc) 00187 {} 00188 00189 /** construct with given shape and copy data from C-style array \a init. 00190 Unless \a layout is <tt>ColumnMajor</tt>, the elements in this array 00191 are assumed to be given in row-major order (the C standard order) and 00192 will automatically be converted to the required column-major format. 00193 Note that the order of the axes is <tt>difference_type(rows, columns)</tt> which 00194 is the opposite of the usual VIGRA convention. 00195 */ 00196 Matrix(const difference_type &shape, const_pointer init, RawArrayMemoryLayout layout = RowMajor, 00197 allocator_type const & alloc = allocator_type()) 00198 : BaseType(shape, alloc) // FIXME: this function initializes the memory twice 00199 { 00200 if(layout == RowMajor) 00201 { 00202 difference_type trans(shape[1], shape[0]); 00203 linalg::transpose(MultiArrayView<2, T>(trans, const_cast<pointer>(init)), *this); 00204 } 00205 else 00206 { 00207 std::copy(init, init + elementCount(), this->data()); 00208 } 00209 } 00210 00211 /** construct with given shape and copy data from C-style array \a init. 00212 Unless \a layout is <tt>ColumnMajor</tt>, the elements in this array 00213 are assumed to be given in row-major order (the C standard order) and 00214 will automatically be converted to the required column-major format. 00215 Note that the order of the axes is <tt>(rows, columns)</tt> which 00216 is the opposite of the usual VIGRA convention. 00217 */ 00218 Matrix(difference_type_1 rows, difference_type_1 columns, const_pointer init, RawArrayMemoryLayout layout = RowMajor, 00219 allocator_type const & alloc = allocator_type()) 00220 : BaseType(difference_type(rows, columns), alloc) // FIXME: this function initializes the memory twice 00221 { 00222 if(layout == RowMajor) 00223 { 00224 difference_type trans(columns, rows); 00225 linalg::transpose(MultiArrayView<2, T>(trans, const_cast<pointer>(init)), *this); 00226 } 00227 else 00228 { 00229 std::copy(init, init + elementCount(), this->data()); 00230 } 00231 } 00232 00233 /** copy constructor. Allocates new memory and 00234 copies tha data. 00235 */ 00236 Matrix(const Matrix &rhs) 00237 : BaseType(rhs) 00238 {} 00239 00240 /** construct from temporary matrix, which looses its data. 00241 00242 This operation is equivalent to 00243 \code 00244 TemporaryMatrix<T> temp = ...; 00245 00246 Matrix<T> m; 00247 m.swap(temp); 00248 \endcode 00249 */ 00250 Matrix(const TemporaryMatrix<T, ALLOC> &rhs) 00251 : BaseType(rhs.allocator()) 00252 { 00253 this->swap(const_cast<TemporaryMatrix<T, ALLOC> &>(rhs)); 00254 } 00255 00256 /** construct from a MultiArrayView. Allocates new memory and 00257 copies tha data. \a rhs is assumed to be in column-major order already. 00258 */ 00259 template<class U, class C> 00260 Matrix(const MultiArrayView<2, U, C> &rhs) 00261 : BaseType(rhs) 00262 {} 00263 00264 /** assignment. 00265 If the size of \a rhs is the same as the matrix's old size, only the data 00266 are copied. Otherwise, new storage is allocated, which invalidates 00267 all objects (array views, iterators) depending on the matrix. 00268 */ 00269 Matrix & operator=(const Matrix &rhs) 00270 { 00271 BaseType::operator=(rhs); // has the correct semantics already 00272 return *this; 00273 } 00274 00275 /** assign a temporary matrix. If the shapes of the two matrices match, 00276 only the data are copied (in order to not invalidate views and iterators 00277 depending on this matrix). Otherwise, the memory is swapped 00278 between the two matrices, so that all depending objects 00279 (array views, iterators) ar invalidated. 00280 */ 00281 Matrix & operator=(const TemporaryMatrix<T, ALLOC> &rhs) 00282 { 00283 if(this->shape() == rhs.shape()) 00284 this->copy(rhs); 00285 else 00286 this->swap(const_cast<TemporaryMatrix<T, ALLOC> &>(rhs)); 00287 return *this; 00288 } 00289 00290 /** assignment from arbitrary 2-dimensional MultiArrayView.<br> 00291 If the size of \a rhs is the same as the matrix's old size, only the data 00292 are copied. Otherwise, new storage is allocated, which invalidates 00293 all objects (array views, iterators) depending on the matrix. 00294 \a rhs is assumed to be in column-major order already. 00295 */ 00296 template <class U, class C> 00297 Matrix & operator=(const MultiArrayView<2, U, C> &rhs) 00298 { 00299 BaseType::operator=(rhs); // has the correct semantics already 00300 return *this; 00301 } 00302 00303 /** init elements with a constant 00304 */ 00305 template <class U> 00306 Matrix & init(const U & init) 00307 { 00308 BaseType::init(init); 00309 return *this; 00310 } 00311 00312 /** reshape to the given shape and initialize with zero. 00313 */ 00314 void reshape(difference_type_1 rows, difference_type_1 columns) 00315 { 00316 BaseType::reshape(difference_type(rows, columns)); 00317 } 00318 00319 /** reshape to the given shape and initialize with \a init. 00320 */ 00321 void reshape(difference_type_1 rows, difference_type_1 columns, const_reference init) 00322 { 00323 BaseType::reshape(difference_type(rows, columns), init); 00324 } 00325 00326 /** reshape to the given shape and initialize with zero. 00327 */ 00328 void reshape(difference_type const & shape) 00329 { 00330 BaseType::reshape(shape); 00331 } 00332 00333 /** reshape to the given shape and initialize with \a init. 00334 */ 00335 void reshape(difference_type const & shape, const_reference init) 00336 { 00337 BaseType::reshape(shape, init); 00338 } 00339 00340 /** Create a matrix view that represents the row vector of row \a d. 00341 */ 00342 view_type rowVector(difference_type_1 d) const 00343 { 00344 return vigra::linalg::rowVector(*this, d); 00345 } 00346 00347 /** Create a matrix view that represents the column vector of column \a d. 00348 */ 00349 view_type columnVector(difference_type_1 d) const 00350 { 00351 return vigra::linalg::columnVector(*this, d); 00352 } 00353 00354 /** number of rows (height) of the matrix. 00355 */ 00356 difference_type_1 rowCount() const 00357 { 00358 return this->m_shape[0]; 00359 } 00360 00361 /** number of columns (width) of the matrix. 00362 */ 00363 difference_type_1 columnCount() const 00364 { 00365 return this->m_shape[1]; 00366 } 00367 00368 /** number of elements (width*height) of the matrix. 00369 */ 00370 difference_type_1 elementCount() const 00371 { 00372 return rowCount()*columnCount(); 00373 } 00374 00375 /** check whether the matrix is symmetric. 00376 */ 00377 bool isSymmetric() const 00378 { 00379 return vigra::linalg::isSymmetric(*this); 00380 } 00381 00382 #ifdef DOXYGEN 00383 // repeat the following functions for documentation. In real code, they are inherited. 00384 00385 /** read/write access to matrix element <tt>(row, column)</tt>. 00386 Note that the order of the argument is the opposite of the usual 00387 VIGRA convention due to column-major matrix order. 00388 */ 00389 value_type & operator()(difference_type_1 row, difference_type_1 column); 00390 00391 /** read access to matrix element <tt>(row, column)</tt>. 00392 Note that the order of the argument is the opposite of the usual 00393 VIGRA convention due to column-major matrix order. 00394 */ 00395 value_type operator()(difference_type_1 row, difference_type_1 column) const; 00396 00397 /** squared Frobenius norm. Sum of squares of the matrix elements. 00398 */ 00399 typename NormTraits<Matrix>::SquaredNormType squaredNorm() const; 00400 00401 /** Frobenius norm. Root of sum of squares of the matrix elements. 00402 */ 00403 typename NormTraits<Matrix>::NormType norm() const; 00404 00405 /** create a transposed view of this matrix. 00406 No data are copied. If you want to transpose this matrix permanently, 00407 you have to assign the transposed view: 00408 00409 \code 00410 a = a.transpose(); 00411 \endcode 00412 */ 00413 MultiArrayView<2, vluae_type, StridedArrayTag> transpose() const; 00414 #endif 00415 00416 /** add \a other to this (sizes must match). 00417 */ 00418 template <class U, class C> 00419 Matrix & operator+=(MultiArrayView<2, U, C> const & other) 00420 { 00421 BaseType::operator+=(other); 00422 return *this; 00423 } 00424 00425 /** subtract \a other from this (sizes must match). 00426 */ 00427 template <class U, class C> 00428 Matrix & operator-=(MultiArrayView<2, U, C> const & other) 00429 { 00430 BaseType::operator-=(other); 00431 return *this; 00432 } 00433 00434 /** multiply \a other element-wise with this matrix (sizes must match). 00435 */ 00436 template <class U, class C> 00437 Matrix & operator*=(MultiArrayView<2, U, C> const & other) 00438 { 00439 BaseType::operator*=(other); 00440 return *this; 00441 } 00442 00443 /** divide this matrix element-wise by \a other (sizes must match). 00444 */ 00445 template <class U, class C> 00446 Matrix & operator/=(MultiArrayView<2, U, C> const & other) 00447 { 00448 BaseType::operator/=(other); 00449 return *this; 00450 } 00451 00452 /** add \a other to each element of this matrix 00453 */ 00454 Matrix & operator+=(T other) 00455 { 00456 BaseType::operator+=(other); 00457 return *this; 00458 } 00459 00460 /** subtraxt \a other from each element of this matrix 00461 */ 00462 Matrix & operator-=(T other) 00463 { 00464 BaseType::operator-=(other); 00465 return *this; 00466 } 00467 00468 /** scalar multiply this with \a other 00469 */ 00470 Matrix & operator*=(T other) 00471 { 00472 BaseType::operator*=(other); 00473 return *this; 00474 } 00475 00476 /** scalar devide this by \a other 00477 */ 00478 Matrix & operator/=(T other) 00479 { 00480 BaseType::operator/=(other); 00481 return *this; 00482 } 00483 }; 00484 00485 // TemporaryMatrix is provided as an optimization: Functions returning a matrix can 00486 // use TemporaryMatrix to make explicit that it was allocated as a temporary data structure. 00487 // Functions receiving a TemporaryMatrix can thus often avoid to allocate new temporary 00488 // memory. 00489 template <class T, class ALLOC = std::allocator<T> > 00490 class TemporaryMatrix 00491 : public Matrix<T, ALLOC> 00492 { 00493 typedef Matrix<T, ALLOC> BaseType; 00494 public: 00495 typedef Matrix<T, ALLOC> matrix_type; 00496 typedef TemporaryMatrix<T, ALLOC> temp_type; 00497 typedef MultiArrayView<2, T, UnstridedArrayTag> view_type; 00498 typedef typename BaseType::value_type value_type; 00499 typedef typename BaseType::pointer pointer; 00500 typedef typename BaseType::const_pointer const_pointer; 00501 typedef typename BaseType::reference reference; 00502 typedef typename BaseType::const_reference const_reference; 00503 typedef typename BaseType::difference_type difference_type; 00504 typedef typename BaseType::difference_type_1 difference_type_1; 00505 typedef ALLOC allocator_type; 00506 00507 TemporaryMatrix(difference_type const & shape) 00508 : BaseType(shape, ALLOC()) 00509 {} 00510 00511 TemporaryMatrix(difference_type const & shape, const_reference init) 00512 : BaseType(shape, init, ALLOC()) 00513 {} 00514 00515 TemporaryMatrix(difference_type_1 rows, difference_type_1 columns) 00516 : BaseType(rows, columns, ALLOC()) 00517 {} 00518 00519 TemporaryMatrix(difference_type_1 rows, difference_type_1 columns, const_reference init) 00520 : BaseType(rows, columns, init, ALLOC()) 00521 {} 00522 00523 template<class U, class C> 00524 TemporaryMatrix(const MultiArrayView<2, U, C> &rhs) 00525 : BaseType(rhs) 00526 {} 00527 00528 TemporaryMatrix(const TemporaryMatrix &rhs) 00529 : BaseType() 00530 { 00531 this->swap(const_cast<TemporaryMatrix &>(rhs)); 00532 } 00533 00534 template <class U> 00535 TemporaryMatrix & init(const U & init) 00536 { 00537 BaseType::init(init); 00538 return *this; 00539 } 00540 00541 template <class U, class C> 00542 TemporaryMatrix & operator+=(MultiArrayView<2, U, C> const & other) 00543 { 00544 BaseType::operator+=(other); 00545 return *this; 00546 } 00547 00548 template <class U, class C> 00549 TemporaryMatrix & operator-=(MultiArrayView<2, U, C> const & other) 00550 { 00551 BaseType::operator-=(other); 00552 return *this; 00553 } 00554 00555 template <class U, class C> 00556 TemporaryMatrix & operator*=(MultiArrayView<2, U, C> const & other) 00557 { 00558 BaseType::operator*=(other); 00559 return *this; 00560 } 00561 00562 template <class U, class C> 00563 TemporaryMatrix & operator/=(MultiArrayView<2, U, C> const & other) 00564 { 00565 BaseType::operator/=(other); 00566 return *this; 00567 } 00568 00569 TemporaryMatrix & operator+=(T other) 00570 { 00571 BaseType::operator+=(other); 00572 return *this; 00573 } 00574 00575 TemporaryMatrix & operator-=(T other) 00576 { 00577 BaseType::operator-=(other); 00578 return *this; 00579 } 00580 00581 TemporaryMatrix & operator*=(T other) 00582 { 00583 BaseType::operator*=(other); 00584 return *this; 00585 } 00586 00587 TemporaryMatrix & operator/=(T other) 00588 { 00589 BaseType::operator/=(other); 00590 return *this; 00591 } 00592 private: 00593 00594 TemporaryMatrix &operator=(const TemporaryMatrix &rhs); // intentionally not implemented 00595 }; 00596 00597 /** \defgroup LinearAlgebraFunctions Matrix Functions 00598 00599 \brief Basic matrix algebra, element-wise mathematical functions, row and columns statistics, data normalization etc. 00600 00601 \ingroup LinearAlgebraModule 00602 */ 00603 //@{ 00604 00605 /** Number of rows of a matrix represented as a <tt>MultiArrayView<2, ...></tt> 00606 00607 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 00608 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00609 Namespaces: vigra and vigra::linalg 00610 */ 00611 template <class T, class C> 00612 inline MultiArrayIndex 00613 rowCount(const MultiArrayView<2, T, C> &x) 00614 { 00615 return x.shape(0); 00616 } 00617 00618 /** Number of columns of a matrix represented as a <tt>MultiArrayView<2, ...></tt> 00619 00620 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 00621 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00622 Namespaces: vigra and vigra::linalg 00623 */ 00624 template <class T, class C> 00625 inline MultiArrayIndex 00626 columnCount(const MultiArrayView<2, T, C> &x) 00627 { 00628 return x.shape(1); 00629 } 00630 00631 /** Create a row vector view for row \a d of the matrix \a m 00632 00633 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 00634 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00635 Namespaces: vigra and vigra::linalg 00636 */ 00637 template <class T, class C> 00638 inline MultiArrayView <2, T, C> 00639 rowVector(MultiArrayView <2, T, C> const & m, MultiArrayIndex d) 00640 { 00641 typedef typename MultiArrayView <2, T, C>::difference_type Shape; 00642 return m.subarray(Shape(d, 0), Shape(d+1, columnCount(m))); 00643 } 00644 00645 00646 /** Create a row vector view of the matrix \a m starting at element \a first and ranging 00647 to column \a end (non-inclusive). 00648 00649 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 00650 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00651 Namespaces: vigra and vigra::linalg 00652 */ 00653 template <class T, class C> 00654 inline MultiArrayView <2, T, C> 00655 rowVector(MultiArrayView <2, T, C> const & m, MultiArrayShape<2>::type first, MultiArrayIndex end) 00656 { 00657 typedef typename MultiArrayView <2, T, C>::difference_type Shape; 00658 return m.subarray(first, Shape(first[0]+1, end)); 00659 } 00660 00661 /** Create a column vector view for column \a d of the matrix \a m 00662 00663 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 00664 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00665 Namespaces: vigra and vigra::linalg 00666 */ 00667 template <class T, class C> 00668 inline MultiArrayView <2, T, C> 00669 columnVector(MultiArrayView<2, T, C> const & m, MultiArrayIndex d) 00670 { 00671 typedef typename MultiArrayView <2, T, C>::difference_type Shape; 00672 return m.subarray(Shape(0, d), Shape(rowCount(m), d+1)); 00673 } 00674 00675 /** Create a column vector view of the matrix \a m starting at element \a first and 00676 ranging to row \a end (non-inclusive). 00677 00678 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 00679 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00680 Namespaces: vigra and vigra::linalg 00681 **/ 00682 template <class T, class C> 00683 inline MultiArrayView <2, T, C> 00684 columnVector(MultiArrayView<2, T, C> const & m, MultiArrayShape<2>::type first, int end) 00685 { 00686 typedef typename MultiArrayView <2, T, C>::difference_type Shape; 00687 return m.subarray(first, Shape(end, first[1]+1)); 00688 } 00689 00690 /** Create a sub vector view of the vector \a m starting at element \a first and 00691 ranging to row \a end (non-inclusive). 00692 00693 Note: This function may only be called when either <tt>rowCount(m) == 1</tt> or 00694 <tt>columnCount(m) == 1</tt>, i.e. when \a m really represents a vector. 00695 Otherwise, a PreconditionViolation exception is raised. 00696 00697 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 00698 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00699 Namespaces: vigra and vigra::linalg 00700 **/ 00701 template <class T, class C> 00702 inline MultiArrayView <2, T, C> 00703 subVector(MultiArrayView<2, T, C> const & m, int first, int end) 00704 { 00705 typedef typename MultiArrayView <2, T, C>::difference_type Shape; 00706 if(columnCount(m) == 1) 00707 return m.subarray(Shape(first, 0), Shape(end, 1)); 00708 vigra_precondition(rowCount(m) == 1, 00709 "linalg::subVector(): Input must be a vector (1xN or Nx1)."); 00710 return m.subarray(Shape(0, first), Shape(1, end)); 00711 } 00712 00713 /** Check whether matrix \a m is symmetric. 00714 00715 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 00716 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00717 Namespaces: vigra and vigra::linalg 00718 */ 00719 template <class T, class C> 00720 bool 00721 isSymmetric(MultiArrayView<2, T, C> const & m) 00722 { 00723 const MultiArrayIndex size = rowCount(m); 00724 if(size != columnCount(m)) 00725 return false; 00726 00727 for(MultiArrayIndex i = 0; i < size; ++i) 00728 for(MultiArrayIndex j = i+1; j < size; ++j) 00729 if(m(j, i) != m(i, j)) 00730 return false; 00731 return true; 00732 } 00733 00734 00735 /** Compute the trace of a square matrix. 00736 00737 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 00738 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00739 Namespaces: vigra and vigra::linalg 00740 */ 00741 template <class T, class C> 00742 typename NumericTraits<T>::Promote 00743 trace(MultiArrayView<2, T, C> const & m) 00744 { 00745 typedef typename NumericTraits<T>::Promote SumType; 00746 00747 const MultiArrayIndex size = rowCount(m); 00748 vigra_precondition(size == columnCount(m), "linalg::trace(): Matrix must be square."); 00749 00750 SumType sum = NumericTraits<SumType>::zero(); 00751 for(MultiArrayIndex i = 0; i < size; ++i) 00752 sum += m(i, i); 00753 return sum; 00754 } 00755 00756 #ifdef DOXYGEN // documentation only -- function is already defined in vigra/multi_array.hxx 00757 00758 /** calculate the squared Frobenius norm of a matrix. 00759 Equal to the sum of squares of the matrix elements. 00760 00761 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> 00762 Namespace: vigra 00763 */ 00764 template <class T, class ALLOC> 00765 typename Matrix<T, ALLLOC>::SquaredNormType 00766 squaredNorm(const Matrix<T, ALLLOC> &a); 00767 00768 /** calculate the Frobenius norm of a matrix. 00769 Equal to the root of the sum of squares of the matrix elements. 00770 00771 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> 00772 Namespace: vigra 00773 */ 00774 template <class T, class ALLOC> 00775 typename Matrix<T, ALLLOC>::NormType 00776 norm(const Matrix<T, ALLLOC> &a); 00777 00778 #endif // DOXYGEN 00779 00780 /** initialize the given square matrix as an identity matrix. 00781 00782 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 00783 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00784 Namespaces: vigra and vigra::linalg 00785 */ 00786 template <class T, class C> 00787 void identityMatrix(MultiArrayView<2, T, C> &r) 00788 { 00789 const MultiArrayIndex rows = rowCount(r); 00790 vigra_precondition(rows == columnCount(r), 00791 "identityMatrix(): Matrix must be square."); 00792 for(MultiArrayIndex i = 0; i < rows; ++i) { 00793 for(MultiArrayIndex j = 0; j < rows; ++j) 00794 r(j, i) = NumericTraits<T>::zero(); 00795 r(i, i) = NumericTraits<T>::one(); 00796 } 00797 } 00798 00799 /** create n identity matrix of the given size. 00800 Usage: 00801 00802 \code 00803 vigra::Matrix<double> m = vigra::identityMatrix<double>(size); 00804 \endcode 00805 00806 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 00807 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00808 Namespaces: vigra and vigra::linalg 00809 */ 00810 template <class T> 00811 TemporaryMatrix<T> identityMatrix(MultiArrayIndex size) 00812 { 00813 TemporaryMatrix<T> ret(size, size, NumericTraits<T>::zero()); 00814 for(MultiArrayIndex i = 0; i < size; ++i) 00815 ret(i, i) = NumericTraits<T>::one(); 00816 return ret; 00817 } 00818 00819 template <class T, class C1, class C2> 00820 void diagonalMatrixImpl(MultiArrayView<1, T, C1> const & v, MultiArrayView<2, T, C2> &r) 00821 { 00822 const MultiArrayIndex size = v.elementCount(); 00823 vigra_precondition(rowCount(r) == size && columnCount(r) == size, 00824 "diagonalMatrix(): result must be a square matrix."); 00825 for(MultiArrayIndex i = 0; i < size; ++i) 00826 r(i, i) = v(i); 00827 } 00828 00829 /** make a diagonal matrix from a vector. 00830 The vector is given as matrix \a v, which must either have a single 00831 row or column. The result is witten into the square matrix \a r. 00832 00833 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 00834 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00835 Namespaces: vigra and vigra::linalg 00836 */ 00837 template <class T, class C1, class C2> 00838 void diagonalMatrix(MultiArrayView<2, T, C1> const & v, MultiArrayView<2, T, C2> &r) 00839 { 00840 vigra_precondition(rowCount(v) == 1 || columnCount(v) == 1, 00841 "diagonalMatrix(): input must be a vector."); 00842 r.init(NumericTraits<T>::zero()); 00843 if(rowCount(v) == 1) 00844 diagonalMatrixImpl(v.bindInner(0), r); 00845 else 00846 diagonalMatrixImpl(v.bindOuter(0), r); 00847 } 00848 00849 /** create a diagonal matrix from a vector. 00850 The vector is given as matrix \a v, which must either have a single 00851 row or column. The result is returned as a temporary matrix. 00852 Usage: 00853 00854 \code 00855 vigra::Matrix<double> v(1, len); 00856 v = ...; 00857 00858 vigra::Matrix<double> m = diagonalMatrix(v); 00859 \endcode 00860 00861 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 00862 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00863 Namespaces: vigra and vigra::linalg 00864 */ 00865 template <class T, class C> 00866 TemporaryMatrix<T> diagonalMatrix(MultiArrayView<2, T, C> const & v) 00867 { 00868 vigra_precondition(rowCount(v) == 1 || columnCount(v) == 1, 00869 "diagonalMatrix(): input must be a vector."); 00870 MultiArrayIndex size = v.elementCount(); 00871 TemporaryMatrix<T> ret(size, size, NumericTraits<T>::zero()); 00872 if(rowCount(v) == 1) 00873 diagonalMatrixImpl(v.bindInner(0), ret); 00874 else 00875 diagonalMatrixImpl(v.bindOuter(0), ret); 00876 return ret; 00877 } 00878 00879 /** transpose matrix \a v. 00880 The result is written into \a r which must have the correct (i.e. 00881 transposed) shape. 00882 00883 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 00884 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00885 Namespaces: vigra and vigra::linalg 00886 */ 00887 template <class T, class C1, class C2> 00888 void transpose(const MultiArrayView<2, T, C1> &v, MultiArrayView<2, T, C2> &r) 00889 { 00890 const MultiArrayIndex rows = rowCount(r); 00891 const MultiArrayIndex cols = columnCount(r); 00892 vigra_precondition(rows == columnCount(v) && cols == rowCount(v), 00893 "transpose(): arrays must have transposed shapes."); 00894 for(MultiArrayIndex i = 0; i < cols; ++i) 00895 for(MultiArrayIndex j = 0; j < rows; ++j) 00896 r(j, i) = v(i, j); 00897 } 00898 00899 /** create the transpose of matrix \a v. 00900 This does not copy any data, but only creates a transposed view 00901 to the original matrix. A copy is only made when the transposed view 00902 is assigned to another matrix. 00903 Usage: 00904 00905 \code 00906 vigra::Matrix<double> v(rows, cols); 00907 v = ...; 00908 00909 vigra::Matrix<double> m = transpose(v); 00910 \endcode 00911 00912 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 00913 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00914 Namespaces: vigra and vigra::linalg 00915 */ 00916 template <class T, class C> 00917 inline MultiArrayView<2, T, StridedArrayTag> 00918 transpose(MultiArrayView<2, T, C> const & v) 00919 { 00920 return v.transpose(); 00921 } 00922 00923 /** Create new matrix by concatenating two matrices \a a and \a b vertically, i.e. on top of each other. 00924 The two matrices must have the same number of columns. 00925 The result is returned as a temporary matrix. 00926 00927 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 00928 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00929 Namespace: vigra::linalg 00930 */ 00931 template <class T, class C1, class C2> 00932 inline TemporaryMatrix<T> 00933 joinVertically(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b) 00934 { 00935 typedef typename TemporaryMatrix<T>::difference_type Shape; 00936 00937 MultiArrayIndex n = columnCount(a); 00938 vigra_precondition(n == columnCount(b), 00939 "joinVertically(): shape mismatch."); 00940 00941 MultiArrayIndex ma = rowCount(a); 00942 MultiArrayIndex mb = rowCount(b); 00943 TemporaryMatrix<T> t(ma + mb, n, T()); 00944 t.subarray(Shape(0,0), Shape(ma, n)) = a; 00945 t.subarray(Shape(ma,0), Shape(ma+mb, n)) = b; 00946 return t; 00947 } 00948 00949 /** Create new matrix by concatenating two matrices \a a and \a b horizontally, i.e. side by side. 00950 The two matrices must have the same number of rows. 00951 The result is returned as a temporary matrix. 00952 00953 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 00954 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00955 Namespace: vigra::linalg 00956 */ 00957 template <class T, class C1, class C2> 00958 inline TemporaryMatrix<T> 00959 joinHorizontally(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b) 00960 { 00961 typedef typename TemporaryMatrix<T>::difference_type Shape; 00962 00963 MultiArrayIndex m = rowCount(a); 00964 vigra_precondition(m == rowCount(b), 00965 "joinHorizontally(): shape mismatch."); 00966 00967 MultiArrayIndex na = columnCount(a); 00968 MultiArrayIndex nb = columnCount(b); 00969 TemporaryMatrix<T> t(m, na + nb, T()); 00970 t.subarray(Shape(0,0), Shape(m, na)) = a; 00971 t.subarray(Shape(0, na), Shape(m, na + nb)) = b; 00972 return t; 00973 } 00974 00975 /** Initialize a matrix with repeated copies of a given matrix. 00976 00977 Matrix \a r will consist of \a verticalCount downward repetitions of \a v, 00978 and \a horizontalCount side-by-side repetitions. When \a v has size <tt>m</tt> by <tt>n</tt>, 00979 \a r must have size <tt>(m*verticalCount)</tt> by <tt>(n*horizontalCount)</tt>. 00980 00981 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 00982 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00983 Namespace: vigra::linalg 00984 */ 00985 template <class T, class C1, class C2> 00986 void repeatMatrix(MultiArrayView<2, T, C1> const & v, MultiArrayView<2, T, C2> &r, 00987 unsigned int verticalCount, unsigned int horizontalCount) 00988 { 00989 typedef typename Matrix<T>::difference_type Shape; 00990 00991 MultiArrayIndex m = rowCount(v), n = columnCount(v); 00992 vigra_precondition(m*verticalCount == rowCount(r) && n*horizontalCount == columnCount(r), 00993 "repeatMatrix(): Shape mismatch."); 00994 00995 for(MultiArrayIndex l=0; l<(MultiArrayIndex)horizontalCount; ++l) 00996 { 00997 for(MultiArrayIndex k=0; k<(MultiArrayIndex)verticalCount; ++k) 00998 { 00999 r.subarray(Shape(k*m, l*n), Shape((k+1)*m, (l+1)*n)) = v; 01000 } 01001 } 01002 } 01003 01004 /** Create a new matrix by repeating a given matrix. 01005 01006 The resulting matrix \a r will consist of \a verticalCount downward repetitions of \a v, 01007 and \a horizontalCount side-by-side repetitions, i.e. it will be of size 01008 <tt>(m*verticalCount)</tt> by <tt>(n*horizontalCount)</tt> when \a v has size <tt>m</tt> by <tt>n</tt>. 01009 The result is returned as a temporary matrix. 01010 01011 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01012 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01013 Namespace: vigra::linalg 01014 */ 01015 template <class T, class C> 01016 TemporaryMatrix<T> 01017 repeatMatrix(MultiArrayView<2, T, C> const & v, unsigned int verticalCount, unsigned int horizontalCount) 01018 { 01019 MultiArrayIndex m = rowCount(v), n = columnCount(v); 01020 TemporaryMatrix<T> ret(verticalCount*m, horizontalCount*n); 01021 repeatMatrix(v, ret, verticalCount, horizontalCount); 01022 return ret; 01023 } 01024 01025 /** add matrices \a a and \a b. 01026 The result is written into \a r. All three matrices must have the same shape. 01027 01028 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01029 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01030 Namespace: vigra::linalg 01031 */ 01032 template <class T, class C1, class C2, class C3> 01033 void add(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b, 01034 MultiArrayView<2, T, C3> &r) 01035 { 01036 const MultiArrayIndex rrows = rowCount(r); 01037 const MultiArrayIndex rcols = columnCount(r); 01038 vigra_precondition(rrows == rowCount(a) && rcols == columnCount(a) && 01039 rrows == rowCount(b) && rcols == columnCount(b), 01040 "add(): Matrix shapes must agree."); 01041 01042 for(MultiArrayIndex i = 0; i < rcols; ++i) { 01043 for(MultiArrayIndex j = 0; j < rrows; ++j) { 01044 r(j, i) = a(j, i) + b(j, i); 01045 } 01046 } 01047 } 01048 01049 /** add matrices \a a and \a b. 01050 The two matrices must have the same shape. 01051 The result is returned as a temporary matrix. 01052 01053 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01054 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01055 Namespace: vigra::linalg 01056 */ 01057 template <class T, class C1, class C2> 01058 inline TemporaryMatrix<T> 01059 operator+(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b) 01060 { 01061 return TemporaryMatrix<T>(a) += b; 01062 } 01063 01064 template <class T, class C> 01065 inline TemporaryMatrix<T> 01066 operator+(const TemporaryMatrix<T> &a, const MultiArrayView<2, T, C> &b) 01067 { 01068 return const_cast<TemporaryMatrix<T> &>(a) += b; 01069 } 01070 01071 template <class T, class C> 01072 inline TemporaryMatrix<T> 01073 operator+(const MultiArrayView<2, T, C> &a, const TemporaryMatrix<T> &b) 01074 { 01075 return const_cast<TemporaryMatrix<T> &>(b) += a; 01076 } 01077 01078 template <class T> 01079 inline TemporaryMatrix<T> 01080 operator+(const TemporaryMatrix<T> &a, const TemporaryMatrix<T> &b) 01081 { 01082 return const_cast<TemporaryMatrix<T> &>(a) += b; 01083 } 01084 01085 /** add scalar \a b to every element of the matrix \a a. 01086 The result is returned as a temporary matrix. 01087 01088 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01089 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01090 Namespace: vigra::linalg 01091 */ 01092 template <class T, class C> 01093 inline TemporaryMatrix<T> 01094 operator+(const MultiArrayView<2, T, C> &a, T b) 01095 { 01096 return TemporaryMatrix<T>(a) += b; 01097 } 01098 01099 template <class T> 01100 inline TemporaryMatrix<T> 01101 operator+(const TemporaryMatrix<T> &a, T b) 01102 { 01103 return const_cast<TemporaryMatrix<T> &>(a) += b; 01104 } 01105 01106 /** add scalar \a a to every element of the matrix \a b. 01107 The result is returned as a temporary matrix. 01108 01109 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01110 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01111 Namespace: vigra::linalg 01112 */ 01113 template <class T, class C> 01114 inline TemporaryMatrix<T> 01115 operator+(T a, const MultiArrayView<2, T, C> &b) 01116 { 01117 return TemporaryMatrix<T>(b) += a; 01118 } 01119 01120 template <class T> 01121 inline TemporaryMatrix<T> 01122 operator+(T a, const TemporaryMatrix<T> &b) 01123 { 01124 return const_cast<TemporaryMatrix<T> &>(b) += a; 01125 } 01126 01127 /** subtract matrix \a b from \a a. 01128 The result is written into \a r. All three matrices must have the same shape. 01129 01130 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01131 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01132 Namespace: vigra::linalg 01133 */ 01134 template <class T, class C1, class C2, class C3> 01135 void sub(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b, 01136 MultiArrayView<2, T, C3> &r) 01137 { 01138 const MultiArrayIndex rrows = rowCount(r); 01139 const MultiArrayIndex rcols = columnCount(r); 01140 vigra_precondition(rrows == rowCount(a) && rcols == columnCount(a) && 01141 rrows == rowCount(b) && rcols == columnCount(b), 01142 "subtract(): Matrix shapes must agree."); 01143 01144 for(MultiArrayIndex i = 0; i < rcols; ++i) { 01145 for(MultiArrayIndex j = 0; j < rrows; ++j) { 01146 r(j, i) = a(j, i) - b(j, i); 01147 } 01148 } 01149 } 01150 01151 /** subtract matrix \a b from \a a. 01152 The two matrices must have the same shape. 01153 The result is returned as a temporary matrix. 01154 01155 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01156 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01157 Namespace: vigra::linalg 01158 */ 01159 template <class T, class C1, class C2> 01160 inline TemporaryMatrix<T> 01161 operator-(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b) 01162 { 01163 return TemporaryMatrix<T>(a) -= b; 01164 } 01165 01166 template <class T, class C> 01167 inline TemporaryMatrix<T> 01168 operator-(const TemporaryMatrix<T> &a, const MultiArrayView<2, T, C> &b) 01169 { 01170 return const_cast<TemporaryMatrix<T> &>(a) -= b; 01171 } 01172 01173 template <class T, class C> 01174 TemporaryMatrix<T> 01175 operator-(const MultiArrayView<2, T, C> &a, const TemporaryMatrix<T> &b) 01176 { 01177 const MultiArrayIndex rows = rowCount(a); 01178 const MultiArrayIndex cols = columnCount(a); 01179 vigra_precondition(rows == b.rowCount() && cols == b.columnCount(), 01180 "Matrix::operator-(): Shape mismatch."); 01181 01182 for(MultiArrayIndex i = 0; i < cols; ++i) 01183 for(MultiArrayIndex j = 0; j < rows; ++j) 01184 const_cast<TemporaryMatrix<T> &>(b)(j, i) = a(j, i) - b(j, i); 01185 return b; 01186 } 01187 01188 template <class T> 01189 inline TemporaryMatrix<T> 01190 operator-(const TemporaryMatrix<T> &a, const TemporaryMatrix<T> &b) 01191 { 01192 return const_cast<TemporaryMatrix<T> &>(a) -= b; 01193 } 01194 01195 /** negate matrix \a a. 01196 The result is returned as a temporary matrix. 01197 01198 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01199 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01200 Namespace: vigra::linalg 01201 */ 01202 template <class T, class C> 01203 inline TemporaryMatrix<T> 01204 operator-(const MultiArrayView<2, T, C> &a) 01205 { 01206 return TemporaryMatrix<T>(a) *= -NumericTraits<T>::one(); 01207 } 01208 01209 template <class T> 01210 inline TemporaryMatrix<T> 01211 operator-(const TemporaryMatrix<T> &a) 01212 { 01213 return const_cast<TemporaryMatrix<T> &>(a) *= -NumericTraits<T>::one(); 01214 } 01215 01216 /** subtract scalar \a b from every element of the matrix \a a. 01217 The result is returned as a temporary matrix. 01218 01219 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01220 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01221 Namespace: vigra::linalg 01222 */ 01223 template <class T, class C> 01224 inline TemporaryMatrix<T> 01225 operator-(const MultiArrayView<2, T, C> &a, T b) 01226 { 01227 return TemporaryMatrix<T>(a) -= b; 01228 } 01229 01230 template <class T> 01231 inline TemporaryMatrix<T> 01232 operator-(const TemporaryMatrix<T> &a, T b) 01233 { 01234 return const_cast<TemporaryMatrix<T> &>(a) -= b; 01235 } 01236 01237 /** subtract every element of the matrix \a b from scalar \a a. 01238 The result is returned as a temporary matrix. 01239 01240 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01241 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01242 Namespace: vigra::linalg 01243 */ 01244 template <class T, class C> 01245 inline TemporaryMatrix<T> 01246 operator-(T a, const MultiArrayView<2, T, C> &b) 01247 { 01248 return TemporaryMatrix<T>(b.shape(), a) -= b; 01249 } 01250 01251 /** calculate the inner product of two matrices representing vectors. 01252 Typically, matrix \a x has a single row, and matrix \a y has 01253 a single column, and the other dimensions match. In addition, this 01254 function handles the cases when either or both of the two inputs are 01255 transposed (e.g. it can compute the dot product of two column vectors). 01256 A <tt>PreconditionViolation</tt> exception is thrown when 01257 the shape conditions are violated. 01258 01259 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01260 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01261 Namespaces: vigra and vigra::linalg 01262 */ 01263 template <class T, class C1, class C2> 01264 typename NormTraits<T>::SquaredNormType 01265 dot(const MultiArrayView<2, T, C1> &x, const MultiArrayView<2, T, C2> &y) 01266 { 01267 typename NormTraits<T>::SquaredNormType ret = 01268 NumericTraits<typename NormTraits<T>::SquaredNormType>::zero(); 01269 if(y.shape(1) == 1) 01270 { 01271 std::ptrdiff_t size = y.shape(0); 01272 if(x.shape(0) == 1 && x.shape(1) == size) // proper scalar product 01273 for(std::ptrdiff_t i = 0; i < size; ++i) 01274 ret += x(0, i) * y(i, 0); 01275 else if(x.shape(1) == 1u && x.shape(0) == size) // two column vectors 01276 for(std::ptrdiff_t i = 0; i < size; ++i) 01277 ret += x(i, 0) * y(i, 0); 01278 else 01279 vigra_precondition(false, "dot(): wrong matrix shapes."); 01280 } 01281 else if(y.shape(0) == 1) 01282 { 01283 std::ptrdiff_t size = y.shape(1); 01284 if(x.shape(0) == 1u && x.shape(1) == size) // two row vectors 01285 for(std::ptrdiff_t i = 0; i < size; ++i) 01286 ret += x(0, i) * y(0, i); 01287 else if(x.shape(1) == 1u && x.shape(0) == size) // column dot row 01288 for(std::ptrdiff_t i = 0; i < size; ++i) 01289 ret += x(i, 0) * y(0, i); 01290 else 01291 vigra_precondition(false, "dot(): wrong matrix shapes."); 01292 } 01293 else 01294 vigra_precondition(false, "dot(): wrong matrix shapes."); 01295 return ret; 01296 } 01297 01298 /** calculate the inner product of two vectors. The vector 01299 lengths must match. 01300 01301 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01302 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01303 Namespaces: vigra and vigra::linalg 01304 */ 01305 template <class T, class C1, class C2> 01306 typename NormTraits<T>::SquaredNormType 01307 dot(const MultiArrayView<1, T, C1> &x, const MultiArrayView<1, T, C2> &y) 01308 { 01309 const MultiArrayIndex n = x.elementCount(); 01310 vigra_precondition(n == y.elementCount(), 01311 "dot(): shape mismatch."); 01312 typename NormTraits<T>::SquaredNormType ret = 01313 NumericTraits<typename NormTraits<T>::SquaredNormType>::zero(); 01314 for(MultiArrayIndex i = 0; i < n; ++i) 01315 ret += x(i) * y(i); 01316 return ret; 01317 } 01318 01319 /** calculate the cross product of two vectors of length 3. 01320 The result is written into \a r. 01321 01322 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01323 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01324 Namespaces: vigra and vigra::linalg 01325 */ 01326 template <class T, class C1, class C2, class C3> 01327 void cross(const MultiArrayView<1, T, C1> &x, const MultiArrayView<1, T, C2> &y, 01328 MultiArrayView<1, T, C3> &r) 01329 { 01330 vigra_precondition(3 == x.elementCount() && 3 == y.elementCount() && 3 == r.elementCount(), 01331 "cross(): vectors must have length 3."); 01332 r(0) = x(1)*y(2) - x(2)*y(1); 01333 r(1) = x(2)*y(0) - x(0)*y(2); 01334 r(2) = x(0)*y(1) - x(1)*y(0); 01335 } 01336 01337 /** calculate the cross product of two matrices representing vectors. 01338 That is, \a x, \a y, and \a r must have a single column of length 3. The result 01339 is written into \a r. 01340 01341 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01342 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01343 Namespaces: vigra and vigra::linalg 01344 */ 01345 template <class T, class C1, class C2, class C3> 01346 void cross(const MultiArrayView<2, T, C1> &x, const MultiArrayView<2, T, C2> &y, 01347 MultiArrayView<2, T, C3> &r) 01348 { 01349 vigra_precondition(3 == rowCount(x) && 3 == rowCount(y) && 3 == rowCount(r) , 01350 "cross(): vectors must have length 3."); 01351 r(0,0) = x(1,0)*y(2,0) - x(2,0)*y(1,0); 01352 r(1,0) = x(2,0)*y(0,0) - x(0,0)*y(2,0); 01353 r(2,0) = x(0,0)*y(1,0) - x(1,0)*y(0,0); 01354 } 01355 01356 /** calculate the cross product of two matrices representing vectors. 01357 That is, \a x, and \a y must have a single column of length 3. The result 01358 is returned as a temporary matrix. 01359 01360 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01361 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01362 Namespaces: vigra and vigra::linalg 01363 */ 01364 template <class T, class C1, class C2> 01365 TemporaryMatrix<T> 01366 cross(const MultiArrayView<2, T, C1> &x, const MultiArrayView<2, T, C2> &y) 01367 { 01368 TemporaryMatrix<T> ret(3, 1); 01369 cross(x, y, ret); 01370 return ret; 01371 } 01372 /** calculate the outer product of two matrices representing vectors. 01373 That is, matrix \a x must have a single column, and matrix \a y must 01374 have a single row, and the other dimensions must match. The result 01375 is written into \a r. 01376 01377 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01378 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01379 Namespaces: vigra and vigra::linalg 01380 */ 01381 template <class T, class C1, class C2, class C3> 01382 void outer(const MultiArrayView<2, T, C1> &x, const MultiArrayView<2, T, C2> &y, 01383 MultiArrayView<2, T, C3> &r) 01384 { 01385 const MultiArrayIndex rows = rowCount(r); 01386 const MultiArrayIndex cols = columnCount(r); 01387 vigra_precondition(rows == rowCount(x) && cols == columnCount(y) && 01388 1 == columnCount(x) && 1 == rowCount(y), 01389 "outer(): shape mismatch."); 01390 for(MultiArrayIndex i = 0; i < cols; ++i) 01391 for(MultiArrayIndex j = 0; j < rows; ++j) 01392 r(j, i) = x(j, 0) * y(0, i); 01393 } 01394 01395 /** calculate the outer product of two matrices representing vectors. 01396 That is, matrix \a x must have a single column, and matrix \a y must 01397 have a single row, and the other dimensions must match. The result 01398 is returned as a temporary matrix. 01399 01400 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01401 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01402 Namespaces: vigra and vigra::linalg 01403 */ 01404 template <class T, class C1, class C2> 01405 TemporaryMatrix<T> 01406 outer(const MultiArrayView<2, T, C1> &x, const MultiArrayView<2, T, C2> &y) 01407 { 01408 const MultiArrayIndex rows = rowCount(x); 01409 const MultiArrayIndex cols = columnCount(y); 01410 vigra_precondition(1 == columnCount(x) && 1 == rowCount(y), 01411 "outer(): shape mismatch."); 01412 TemporaryMatrix<T> ret(rows, cols); 01413 outer(x, y, ret); 01414 return ret; 01415 } 01416 01417 /** calculate the outer product of a matrix (representing a vector) with itself. 01418 The result is returned as a temporary matrix. 01419 01420 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01421 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01422 Namespaces: vigra and vigra::linalg 01423 */ 01424 template <class T, class C> 01425 TemporaryMatrix<T> 01426 outer(const MultiArrayView<2, T, C> &x) 01427 { 01428 const MultiArrayIndex rows = rowCount(x); 01429 const MultiArrayIndex cols = columnCount(x); 01430 vigra_precondition(rows == 1 || cols == 1, 01431 "outer(): matrix does not represent a vector."); 01432 const MultiArrayIndex size = std::max(rows, cols); 01433 TemporaryMatrix<T> ret(size, size); 01434 01435 if(rows == 1) 01436 { 01437 for(MultiArrayIndex i = 0; i < size; ++i) 01438 for(MultiArrayIndex j = 0; j < size; ++j) 01439 ret(j, i) = x(0, j) * x(0, i); 01440 } 01441 else 01442 { 01443 for(MultiArrayIndex i = 0; i < size; ++i) 01444 for(MultiArrayIndex j = 0; j < size; ++j) 01445 ret(j, i) = x(j, 0) * x(i, 0); 01446 } 01447 return ret; 01448 } 01449 01450 template <class T> 01451 class PointWise 01452 { 01453 public: 01454 T const & t; 01455 01456 PointWise(T const & it) 01457 : t(it) 01458 {} 01459 }; 01460 01461 template <class T> 01462 PointWise<T> pointWise(T const & t) 01463 { 01464 return PointWise<T>(t); 01465 } 01466 01467 01468 /** multiply matrix \a a with scalar \a b. 01469 The result is written into \a r. \a a and \a r must have the same shape. 01470 01471 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01472 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01473 Namespace: vigra::linalg 01474 */ 01475 template <class T, class C1, class C2> 01476 void smul(const MultiArrayView<2, T, C1> &a, T b, MultiArrayView<2, T, C2> &r) 01477 { 01478 const MultiArrayIndex rows = rowCount(a); 01479 const MultiArrayIndex cols = columnCount(a); 01480 vigra_precondition(rows == rowCount(r) && cols == columnCount(r), 01481 "smul(): Matrix sizes must agree."); 01482 01483 for(MultiArrayIndex i = 0; i < cols; ++i) 01484 for(MultiArrayIndex j = 0; j < rows; ++j) 01485 r(j, i) = a(j, i) * b; 01486 } 01487 01488 /** multiply scalar \a a with matrix \a b. 01489 The result is written into \a r. \a b and \a r must have the same shape. 01490 01491 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01492 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01493 Namespace: vigra::linalg 01494 */ 01495 template <class T, class C2, class C3> 01496 void smul(T a, const MultiArrayView<2, T, C2> &b, MultiArrayView<2, T, C3> &r) 01497 { 01498 smul(b, a, r); 01499 } 01500 01501 /** perform matrix multiplication of matrices \a a and \a b. 01502 The result is written into \a r. The three matrices must have matching shapes. 01503 01504 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01505 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01506 Namespace: vigra::linalg 01507 */ 01508 template <class T, class C1, class C2, class C3> 01509 void mmul(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b, 01510 MultiArrayView<2, T, C3> &r) 01511 { 01512 const MultiArrayIndex rrows = rowCount(r); 01513 const MultiArrayIndex rcols = columnCount(r); 01514 const MultiArrayIndex acols = columnCount(a); 01515 vigra_precondition(rrows == rowCount(a) && rcols == columnCount(b) && acols == rowCount(b), 01516 "mmul(): Matrix shapes must agree."); 01517 01518 // order of loops ensures that inner loop goes down columns 01519 for(MultiArrayIndex i = 0; i < rcols; ++i) 01520 { 01521 for(MultiArrayIndex j = 0; j < rrows; ++j) 01522 r(j, i) = a(j, 0) * b(0, i); 01523 for(MultiArrayIndex k = 1; k < acols; ++k) 01524 for(MultiArrayIndex j = 0; j < rrows; ++j) 01525 r(j, i) += a(j, k) * b(k, i); 01526 } 01527 } 01528 01529 /** perform matrix multiplication of matrices \a a and \a b. 01530 \a a and \a b must have matching shapes. 01531 The result is returned as a temporary matrix. 01532 01533 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01534 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01535 Namespace: vigra::linalg 01536 */ 01537 template <class T, class C1, class C2> 01538 inline TemporaryMatrix<T> 01539 mmul(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b) 01540 { 01541 TemporaryMatrix<T> ret(rowCount(a), columnCount(b)); 01542 mmul(a, b, ret); 01543 return ret; 01544 } 01545 01546 /** multiply two matrices \a a and \a b pointwise. 01547 The result is written into \a r. All three matrices must have the same shape. 01548 01549 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01550 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01551 Namespace: vigra::linalg 01552 */ 01553 template <class T, class C1, class C2, class C3> 01554 void pmul(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b, 01555 MultiArrayView<2, T, C3> &r) 01556 { 01557 const MultiArrayIndex rrows = rowCount(r); 01558 const MultiArrayIndex rcols = columnCount(r); 01559 vigra_precondition(rrows == rowCount(a) && rcols == columnCount(a) && 01560 rrows == rowCount(b) && rcols == columnCount(b), 01561 "pmul(): Matrix shapes must agree."); 01562 01563 for(MultiArrayIndex i = 0; i < rcols; ++i) { 01564 for(MultiArrayIndex j = 0; j < rrows; ++j) { 01565 r(j, i) = a(j, i) * b(j, i); 01566 } 01567 } 01568 } 01569 01570 /** multiply matrices \a a and \a b pointwise. 01571 \a a and \a b must have matching shapes. 01572 The result is returned as a temporary matrix. 01573 01574 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01575 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01576 Namespace: vigra::linalg 01577 */ 01578 template <class T, class C1, class C2> 01579 inline TemporaryMatrix<T> 01580 pmul(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b) 01581 { 01582 TemporaryMatrix<T> ret(a.shape()); 01583 pmul(a, b, ret); 01584 return ret; 01585 } 01586 01587 /** multiply matrices \a a and \a b pointwise. 01588 \a a and \a b must have matching shapes. 01589 The result is returned as a temporary matrix. 01590 01591 Usage: 01592 01593 \code 01594 Matrix<double> a(m,n), b(m,n); 01595 01596 Matrix<double> c = a * pointWise(b); 01597 // is equivalent to 01598 // Matrix<double> c = pmul(a, b); 01599 \endcode 01600 01601 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01602 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01603 Namespace: vigra::linalg 01604 */ 01605 template <class T, class C, class U> 01606 inline TemporaryMatrix<T> 01607 operator*(const MultiArrayView<2, T, C> &a, PointWise<U> b) 01608 { 01609 return pmul(a, b.t); 01610 } 01611 01612 /** multiply matrix \a a with scalar \a b. 01613 The result is returned as a temporary matrix. 01614 01615 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01616 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01617 Namespace: vigra::linalg 01618 */ 01619 template <class T, class C> 01620 inline TemporaryMatrix<T> 01621 operator*(const MultiArrayView<2, T, C> &a, T b) 01622 { 01623 return TemporaryMatrix<T>(a) *= b; 01624 } 01625 01626 template <class T> 01627 inline TemporaryMatrix<T> 01628 operator*(const TemporaryMatrix<T> &a, T b) 01629 { 01630 return const_cast<TemporaryMatrix<T> &>(a) *= b; 01631 } 01632 01633 /** multiply scalar \a a with matrix \a b. 01634 The result is returned as a temporary matrix. 01635 01636 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01637 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01638 Namespace: vigra::linalg 01639 */ 01640 template <class T, class C> 01641 inline TemporaryMatrix<T> 01642 operator*(T a, const MultiArrayView<2, T, C> &b) 01643 { 01644 return TemporaryMatrix<T>(b) *= a; 01645 } 01646 01647 template <class T> 01648 inline TemporaryMatrix<T> 01649 operator*(T a, const TemporaryMatrix<T> &b) 01650 { 01651 return const_cast<TemporaryMatrix<T> &>(b) *= a; 01652 } 01653 01654 /** multiply matrix \a a with TinyVector \a b. 01655 \a a must be of size <tt>N x N</tt>. Vector \a b and the result 01656 vector are interpreted as column vectors. 01657 01658 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01659 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01660 Namespace: vigra::linalg 01661 */ 01662 template <class T, class A, int N, class DATA, class DERIVED> 01663 TinyVector<T, N> 01664 operator*(const Matrix<T, A> &a, const TinyVectorBase<T, N, DATA, DERIVED> &b) 01665 { 01666 vigra_precondition(N == rowCount(a) && N == columnCount(a), 01667 "operator*(Matrix, TinyVector): Shape mismatch."); 01668 01669 TinyVector<T, N> res = TinyVectorView<T, N>(&a(0,0)) * b[0]; 01670 for(MultiArrayIndex i = 1; i < N; ++i) 01671 res += TinyVectorView<T, N>(&a(0,i)) * b[i]; 01672 return res; 01673 } 01674 01675 /** multiply TinyVector \a a with matrix \a b. 01676 \a b must be of size <tt>N x N</tt>. Vector \a a and the result 01677 vector are interpreted as row vectors. 01678 01679 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01680 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01681 Namespace: vigra::linalg 01682 */ 01683 template <class T, int N, class DATA, class DERIVED, class A> 01684 TinyVector<T, N> 01685 operator*(const TinyVectorBase<T, N, DATA, DERIVED> &a, const Matrix<T, A> &b) 01686 { 01687 vigra_precondition(N == rowCount(b) && N == columnCount(b), 01688 "operator*(TinyVector, Matrix): Shape mismatch."); 01689 01690 TinyVector<T, N> res; 01691 for(MultiArrayIndex i = 0; i < N; ++i) 01692 res[i] = dot(a, TinyVectorView<T, N>(&b(0,i))); 01693 return res; 01694 } 01695 01696 /** perform matrix multiplication of matrices \a a and \a b. 01697 \a a and \a b must have matching shapes. 01698 The result is returned as a temporary matrix. 01699 01700 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01701 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01702 Namespace: vigra::linalg 01703 */ 01704 template <class T, class C1, class C2> 01705 inline TemporaryMatrix<T> 01706 operator*(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b) 01707 { 01708 TemporaryMatrix<T> ret(rowCount(a), columnCount(b)); 01709 mmul(a, b, ret); 01710 return ret; 01711 } 01712 01713 /** divide matrix \a a by scalar \a b. 01714 The result is written into \a r. \a a and \a r must have the same shape. 01715 01716 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01717 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01718 Namespace: vigra::linalg 01719 */ 01720 template <class T, class C1, class C2> 01721 void sdiv(const MultiArrayView<2, T, C1> &a, T b, MultiArrayView<2, T, C2> &r) 01722 { 01723 const MultiArrayIndex rows = rowCount(a); 01724 const MultiArrayIndex cols = columnCount(a); 01725 vigra_precondition(rows == rowCount(r) && cols == columnCount(r), 01726 "sdiv(): Matrix sizes must agree."); 01727 01728 for(MultiArrayIndex i = 0; i < cols; ++i) 01729 for(MultiArrayIndex j = 0; j < rows; ++j) 01730 r(j, i) = a(j, i) / b; 01731 } 01732 01733 /** divide two matrices \a a and \a b pointwise. 01734 The result is written into \a r. All three matrices must have the same shape. 01735 01736 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01737 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01738 Namespace: vigra::linalg 01739 */ 01740 template <class T, class C1, class C2, class C3> 01741 void pdiv(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b, 01742 MultiArrayView<2, T, C3> &r) 01743 { 01744 const MultiArrayIndex rrows = rowCount(r); 01745 const MultiArrayIndex rcols = columnCount(r); 01746 vigra_precondition(rrows == rowCount(a) && rcols == columnCount(a) && 01747 rrows == rowCount(b) && rcols == columnCount(b), 01748 "pdiv(): Matrix shapes must agree."); 01749 01750 for(MultiArrayIndex i = 0; i < rcols; ++i) { 01751 for(MultiArrayIndex j = 0; j < rrows; ++j) { 01752 r(j, i) = a(j, i) / b(j, i); 01753 } 01754 } 01755 } 01756 01757 /** divide matrices \a a and \a b pointwise. 01758 \a a and \a b must have matching shapes. 01759 The result is returned as a temporary matrix. 01760 01761 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01762 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01763 Namespace: vigra::linalg 01764 */ 01765 template <class T, class C1, class C2> 01766 inline TemporaryMatrix<T> 01767 pdiv(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b) 01768 { 01769 TemporaryMatrix<T> ret(a.shape()); 01770 pdiv(a, b, ret); 01771 return ret; 01772 } 01773 01774 /** divide matrices \a a and \a b pointwise. 01775 \a a and \a b must have matching shapes. 01776 The result is returned as a temporary matrix. 01777 01778 Usage: 01779 01780 \code 01781 Matrix<double> a(m,n), b(m,n); 01782 01783 Matrix<double> c = a / pointWise(b); 01784 // is equivalent to 01785 // Matrix<double> c = pdiv(a, b); 01786 \endcode 01787 01788 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01789 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01790 Namespace: vigra::linalg 01791 */ 01792 template <class T, class C, class U> 01793 inline TemporaryMatrix<T> 01794 operator/(const MultiArrayView<2, T, C> &a, PointWise<U> b) 01795 { 01796 return pdiv(a, b.t); 01797 } 01798 01799 /** divide matrix \a a by scalar \a b. 01800 The result is returned as a temporary matrix. 01801 01802 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01803 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01804 Namespace: vigra::linalg 01805 */ 01806 template <class T, class C> 01807 inline TemporaryMatrix<T> 01808 operator/(const MultiArrayView<2, T, C> &a, T b) 01809 { 01810 return TemporaryMatrix<T>(a) /= b; 01811 } 01812 01813 template <class T> 01814 inline TemporaryMatrix<T> 01815 operator/(const TemporaryMatrix<T> &a, T b) 01816 { 01817 return const_cast<TemporaryMatrix<T> &>(a) /= b; 01818 } 01819 01820 /** Create a matrix whose elements are the quotients between scalar \a a and 01821 matrix \a b. The result is returned as a temporary matrix. 01822 01823 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 01824 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01825 Namespace: vigra::linalg 01826 */ 01827 template <class T, class C> 01828 inline TemporaryMatrix<T> 01829 operator/(T a, const MultiArrayView<2, T, C> &b) 01830 { 01831 return TemporaryMatrix<T>(b.shape(), a) / pointWise(b); 01832 } 01833 01834 using vigra::argMin; 01835 using vigra::argMinIf; 01836 using vigra::argMax; 01837 using vigra::argMaxIf; 01838 01839 /*! Find the index of the minimum element in a matrix. 01840 01841 The function returns the index in column-major scan-order sense, 01842 i.e. according to the order used by <tt>MultiArrayView::operator[]</tt>. 01843 If the matrix is actually a vector, this is just the row or columns index. 01844 In case of a truely 2-dimensional matrix, the index can be converted to an 01845 index pair by calling <tt>MultiArrayView::scanOrderIndexToCoordinate()</tt> 01846 01847 <b>Required Interface:</b> 01848 01849 \code 01850 bool f = a[0] < NumericTraits<T>::max(); 01851 \endcode 01852 01853 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>><br> 01854 Namespace: vigra 01855 */ 01856 template <class T, class C> 01857 int argMin(MultiArrayView<2, T, C> const & a) 01858 { 01859 T vopt = NumericTraits<T>::max(); 01860 int best = -1; 01861 for(int k=0; k < a.size(); ++k) 01862 { 01863 if(a[k] < vopt) 01864 { 01865 vopt = a[k]; 01866 best = k; 01867 } 01868 } 01869 return best; 01870 } 01871 01872 /*! Find the index of the maximum element in a matrix. 01873 01874 The function returns the index in column-major scan-order sense, 01875 i.e. according to the order used by <tt>MultiArrayView::operator[]</tt>. 01876 If the matrix is actually a vector, this is just the row or columns index. 01877 In case of a truely 2-dimensional matrix, the index can be converted to an 01878 index pair by calling <tt>MultiArrayView::scanOrderIndexToCoordinate()</tt> 01879 01880 <b>Required Interface:</b> 01881 01882 \code 01883 bool f = NumericTraits<T>::min() < a[0]; 01884 \endcode 01885 01886 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>><br> 01887 Namespace: vigra 01888 */ 01889 template <class T, class C> 01890 int argMax(MultiArrayView<2, T, C> const & a) 01891 { 01892 T vopt = NumericTraits<T>::min(); 01893 int best = -1; 01894 for(int k=0; k < a.size(); ++k) 01895 { 01896 if(vopt < a[k]) 01897 { 01898 vopt = a[k]; 01899 best = k; 01900 } 01901 } 01902 return best; 01903 } 01904 01905 /*! Find the index of the minimum element in a matrix subject to a condition. 01906 01907 The function returns <tt>-1</tt> if no element conforms to \a condition. 01908 Otherwise, the index of the maximum element is returned in column-major scan-order sense, 01909 i.e. according to the order used by <tt>MultiArrayView::operator[]</tt>. 01910 If the matrix is actually a vector, this is just the row or columns index. 01911 In case of a truely 2-dimensional matrix, the index can be converted to an 01912 index pair by calling <tt>MultiArrayView::scanOrderIndexToCoordinate()</tt> 01913 01914 <b>Required Interface:</b> 01915 01916 \code 01917 bool c = condition(a[0]); 01918 bool f = a[0] < NumericTraits<T>::max(); 01919 \endcode 01920 01921 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>><br> 01922 Namespace: vigra 01923 */ 01924 template <class T, class C, class UnaryFunctor> 01925 int argMinIf(MultiArrayView<2, T, C> const & a, UnaryFunctor condition) 01926 { 01927 T vopt = NumericTraits<T>::max(); 01928 int best = -1; 01929 for(int k=0; k < a.size(); ++k) 01930 { 01931 if(condition(a[k]) && a[k] < vopt) 01932 { 01933 vopt = a[k]; 01934 best = k; 01935 } 01936 } 01937 return best; 01938 } 01939 01940 /*! Find the index of the maximum element in a matrix subject to a condition. 01941 01942 The function returns <tt>-1</tt> if no element conforms to \a condition. 01943 Otherwise, the index of the maximum element is returned in column-major scan-order sense, 01944 i.e. according to the order used by <tt>MultiArrayView::operator[]</tt>. 01945 If the matrix is actually a vector, this is just the row or columns index. 01946 In case of a truely 2-dimensional matrix, the index can be converted to an 01947 index pair by calling <tt>MultiArrayView::scanOrderIndexToCoordinate()</tt> 01948 01949 <b>Required Interface:</b> 01950 01951 \code 01952 bool c = condition(a[0]); 01953 bool f = NumericTraits<T>::min() < a[0]; 01954 \endcode 01955 01956 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>><br> 01957 Namespace: vigra 01958 */ 01959 template <class T, class C, class UnaryFunctor> 01960 int argMaxIf(MultiArrayView<2, T, C> const & a, UnaryFunctor condition) 01961 { 01962 T vopt = NumericTraits<T>::min(); 01963 int best = -1; 01964 for(int k=0; k < a.size(); ++k) 01965 { 01966 if(condition(a[k]) && vopt < a[k]) 01967 { 01968 vopt = a[k]; 01969 best = k; 01970 } 01971 } 01972 return best; 01973 } 01974 01975 /** Matrix point-wise power. 01976 */ 01977 template <class T, class C> 01978 linalg::TemporaryMatrix<T> pow(MultiArrayView<2, T, C> const & v, T exponent) 01979 { 01980 linalg::TemporaryMatrix<T> t(v.shape()); 01981 MultiArrayIndex m = rowCount(v), n = columnCount(v); 01982 01983 for(MultiArrayIndex i = 0; i < n; ++i) 01984 for(MultiArrayIndex j = 0; j < m; ++j) 01985 t(j, i) = vigra::pow(v(j, i), exponent); 01986 return t; 01987 } 01988 01989 template <class T> 01990 linalg::TemporaryMatrix<T> pow(linalg::TemporaryMatrix<T> const & v, T exponent) 01991 { 01992 linalg::TemporaryMatrix<T> & t = const_cast<linalg::TemporaryMatrix<T> &>(v); 01993 MultiArrayIndex m = rowCount(t), n = columnCount(t); 01994 01995 for(MultiArrayIndex i = 0; i < n; ++i) 01996 for(MultiArrayIndex j = 0; j < m; ++j) 01997 t(j, i) = vigra::pow(t(j, i), exponent); 01998 return t; 01999 } 02000 02001 template <class T, class C> 02002 linalg::TemporaryMatrix<T> pow(MultiArrayView<2, T, C> const & v, int exponent) 02003 { 02004 linalg::TemporaryMatrix<T> t(v.shape()); 02005 MultiArrayIndex m = rowCount(v), n = columnCount(v); 02006 02007 for(MultiArrayIndex i = 0; i < n; ++i) 02008 for(MultiArrayIndex j = 0; j < m; ++j) 02009 t(j, i) = vigra::pow(v(j, i), exponent); 02010 return t; 02011 } 02012 02013 template <class T> 02014 linalg::TemporaryMatrix<T> pow(linalg::TemporaryMatrix<T> const & v, int exponent) 02015 { 02016 linalg::TemporaryMatrix<T> & t = const_cast<linalg::TemporaryMatrix<T> &>(v); 02017 MultiArrayIndex m = rowCount(t), n = columnCount(t); 02018 02019 for(MultiArrayIndex i = 0; i < n; ++i) 02020 for(MultiArrayIndex j = 0; j < m; ++j) 02021 t(j, i) = vigra::pow(t(j, i), exponent); 02022 return t; 02023 } 02024 02025 template <class C> 02026 linalg::TemporaryMatrix<int> pow(MultiArrayView<2, int, C> const & v, int exponent) 02027 { 02028 linalg::TemporaryMatrix<int> t(v.shape()); 02029 MultiArrayIndex m = rowCount(v), n = columnCount(v); 02030 02031 for(MultiArrayIndex i = 0; i < n; ++i) 02032 for(MultiArrayIndex j = 0; j < m; ++j) 02033 t(j, i) = (int)vigra::pow((double)v(j, i), exponent); 02034 return t; 02035 } 02036 02037 inline 02038 linalg::TemporaryMatrix<int> pow(linalg::TemporaryMatrix<int> const & v, int exponent) 02039 { 02040 linalg::TemporaryMatrix<int> & t = const_cast<linalg::TemporaryMatrix<int> &>(v); 02041 MultiArrayIndex m = rowCount(t), n = columnCount(t); 02042 02043 for(MultiArrayIndex i = 0; i < n; ++i) 02044 for(MultiArrayIndex j = 0; j < m; ++j) 02045 t(j, i) = (int)vigra::pow((double)t(j, i), exponent); 02046 return t; 02047 } 02048 02049 /** Matrix point-wise sqrt. */ 02050 template <class T, class C> 02051 linalg::TemporaryMatrix<T> sqrt(MultiArrayView<2, T, C> const & v); 02052 /** Matrix point-wise exp. */ 02053 template <class T, class C> 02054 linalg::TemporaryMatrix<T> exp(MultiArrayView<2, T, C> const & v); 02055 /** Matrix point-wise log. */ 02056 template <class T, class C> 02057 linalg::TemporaryMatrix<T> log(MultiArrayView<2, T, C> const & v); 02058 /** Matrix point-wise log10. */ 02059 template <class T, class C> 02060 linalg::TemporaryMatrix<T> log10(MultiArrayView<2, T, C> const & v); 02061 /** Matrix point-wise sin. */ 02062 template <class T, class C> 02063 linalg::TemporaryMatrix<T> sin(MultiArrayView<2, T, C> const & v); 02064 /** Matrix point-wise asin. */ 02065 template <class T, class C> 02066 linalg::TemporaryMatrix<T> asin(MultiArrayView<2, T, C> const & v); 02067 /** Matrix point-wise cos. */ 02068 template <class T, class C> 02069 linalg::TemporaryMatrix<T> cos(MultiArrayView<2, T, C> const & v); 02070 /** Matrix point-wise acos. */ 02071 template <class T, class C> 02072 linalg::TemporaryMatrix<T> acos(MultiArrayView<2, T, C> const & v); 02073 /** Matrix point-wise tan. */ 02074 template <class T, class C> 02075 linalg::TemporaryMatrix<T> tan(MultiArrayView<2, T, C> const & v); 02076 /** Matrix point-wise atan. */ 02077 template <class T, class C> 02078 linalg::TemporaryMatrix<T> atan(MultiArrayView<2, T, C> const & v); 02079 /** Matrix point-wise round. */ 02080 template <class T, class C> 02081 linalg::TemporaryMatrix<T> round(MultiArrayView<2, T, C> const & v); 02082 /** Matrix point-wise floor. */ 02083 template <class T, class C> 02084 linalg::TemporaryMatrix<T> floor(MultiArrayView<2, T, C> const & v); 02085 /** Matrix point-wise ceil. */ 02086 template <class T, class C> 02087 linalg::TemporaryMatrix<T> ceil(MultiArrayView<2, T, C> const & v); 02088 /** Matrix point-wise abs. */ 02089 template <class T, class C> 02090 linalg::TemporaryMatrix<T> abs(MultiArrayView<2, T, C> const & v); 02091 /** Matrix point-wise square. */ 02092 template <class T, class C> 02093 linalg::TemporaryMatrix<T> sq(MultiArrayView<2, T, C> const & v); 02094 /** Matrix point-wise sign. */ 02095 template <class T, class C> 02096 linalg::TemporaryMatrix<T> sign(MultiArrayView<2, T, C> const & v); 02097 02098 #define VIGRA_MATRIX_UNARY_FUNCTION(FUNCTION, NAMESPACE) \ 02099 using NAMESPACE::FUNCTION; \ 02100 template <class T, class C> \ 02101 linalg::TemporaryMatrix<T> FUNCTION(MultiArrayView<2, T, C> const & v) \ 02102 { \ 02103 linalg::TemporaryMatrix<T> t(v.shape()); \ 02104 MultiArrayIndex m = rowCount(v), n = columnCount(v); \ 02105 \ 02106 for(MultiArrayIndex i = 0; i < n; ++i) \ 02107 for(MultiArrayIndex j = 0; j < m; ++j) \ 02108 t(j, i) = NAMESPACE::FUNCTION(v(j, i)); \ 02109 return t; \ 02110 } \ 02111 \ 02112 template <class T> \ 02113 linalg::TemporaryMatrix<T> FUNCTION(linalg::Matrix<T> const & v) \ 02114 { \ 02115 linalg::TemporaryMatrix<T> t(v.shape()); \ 02116 MultiArrayIndex m = rowCount(v), n = columnCount(v); \ 02117 \ 02118 for(MultiArrayIndex i = 0; i < n; ++i) \ 02119 for(MultiArrayIndex j = 0; j < m; ++j) \ 02120 t(j, i) = NAMESPACE::FUNCTION(v(j, i)); \ 02121 return t; \ 02122 } \ 02123 \ 02124 template <class T> \ 02125 linalg::TemporaryMatrix<T> FUNCTION(linalg::TemporaryMatrix<T> const & v) \ 02126 { \ 02127 linalg::TemporaryMatrix<T> & t = const_cast<linalg::TemporaryMatrix<T> &>(v); \ 02128 MultiArrayIndex m = rowCount(t), n = columnCount(t); \ 02129 \ 02130 for(MultiArrayIndex i = 0; i < n; ++i) \ 02131 for(MultiArrayIndex j = 0; j < m; ++j) \ 02132 t(j, i) = NAMESPACE::FUNCTION(t(j, i)); \ 02133 return v; \ 02134 }\ 02135 }\ 02136 using linalg::FUNCTION;\ 02137 namespace linalg { 02138 02139 VIGRA_MATRIX_UNARY_FUNCTION(sqrt, std) 02140 VIGRA_MATRIX_UNARY_FUNCTION(exp, std) 02141 VIGRA_MATRIX_UNARY_FUNCTION(log, std) 02142 VIGRA_MATRIX_UNARY_FUNCTION(log10, std) 02143 VIGRA_MATRIX_UNARY_FUNCTION(sin, std) 02144 VIGRA_MATRIX_UNARY_FUNCTION(asin, std) 02145 VIGRA_MATRIX_UNARY_FUNCTION(cos, std) 02146 VIGRA_MATRIX_UNARY_FUNCTION(acos, std) 02147 VIGRA_MATRIX_UNARY_FUNCTION(tan, std) 02148 VIGRA_MATRIX_UNARY_FUNCTION(atan, std) 02149 VIGRA_MATRIX_UNARY_FUNCTION(round, vigra) 02150 VIGRA_MATRIX_UNARY_FUNCTION(floor, vigra) 02151 VIGRA_MATRIX_UNARY_FUNCTION(ceil, vigra) 02152 VIGRA_MATRIX_UNARY_FUNCTION(abs, vigra) 02153 VIGRA_MATRIX_UNARY_FUNCTION(sq, vigra) 02154 VIGRA_MATRIX_UNARY_FUNCTION(sign, vigra) 02155 02156 #undef VIGRA_MATRIX_UNARY_FUNCTION 02157 02158 //@} 02159 02160 } // namespace linalg 02161 02162 using linalg::RowMajor; 02163 using linalg::ColumnMajor; 02164 using linalg::Matrix; 02165 using linalg::identityMatrix; 02166 using linalg::diagonalMatrix; 02167 using linalg::transpose; 02168 using linalg::pointWise; 02169 using linalg::trace; 02170 using linalg::dot; 02171 using linalg::cross; 02172 using linalg::outer; 02173 using linalg::rowCount; 02174 using linalg::columnCount; 02175 using linalg::rowVector; 02176 using linalg::columnVector; 02177 using linalg::subVector; 02178 using linalg::isSymmetric; 02179 using linalg::joinHorizontally; 02180 using linalg::joinVertically; 02181 using linalg::argMin; 02182 using linalg::argMinIf; 02183 using linalg::argMax; 02184 using linalg::argMaxIf; 02185 02186 /********************************************************/ 02187 /* */ 02188 /* NormTraits */ 02189 /* */ 02190 /********************************************************/ 02191 02192 template <class T, class ALLOC> 02193 struct NormTraits<Matrix<T, ALLOC> > 02194 : public NormTraits<MultiArray<2, T, ALLOC> > 02195 { 02196 typedef NormTraits<MultiArray<2, T, ALLOC> > BaseType; 02197 typedef Matrix<T, ALLOC> Type; 02198 typedef typename BaseType::SquaredNormType SquaredNormType; 02199 typedef typename BaseType::NormType NormType; 02200 }; 02201 02202 template <class T, class ALLOC> 02203 struct NormTraits<linalg::TemporaryMatrix<T, ALLOC> > 02204 : public NormTraits<Matrix<T, ALLOC> > 02205 { 02206 typedef NormTraits<Matrix<T, ALLOC> > BaseType; 02207 typedef linalg::TemporaryMatrix<T, ALLOC> Type; 02208 typedef typename BaseType::SquaredNormType SquaredNormType; 02209 typedef typename BaseType::NormType NormType; 02210 }; 02211 02212 } // namespace vigra 02213 02214 namespace std { 02215 02216 /** \addtogroup LinearAlgebraFunctions 02217 */ 02218 //@{ 02219 02220 /** print a matrix \a m to the stream \a s. 02221 02222 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 02223 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 02224 Namespace: std 02225 */ 02226 template <class T, class C> 02227 ostream & 02228 operator<<(ostream & s, const vigra::MultiArrayView<2, T, C> &m) 02229 { 02230 const vigra::MultiArrayIndex rows = vigra::linalg::rowCount(m); 02231 const vigra::MultiArrayIndex cols = vigra::linalg::columnCount(m); 02232 ios::fmtflags flags = s.setf(ios::right | ios::fixed, ios::adjustfield | ios::floatfield); 02233 for(vigra::MultiArrayIndex j = 0; j < rows; ++j) 02234 { 02235 for(vigra::MultiArrayIndex i = 0; i < cols; ++i) 02236 { 02237 s << m(j, i) << " "; 02238 } 02239 s << endl; 02240 } 02241 s.setf(flags); 02242 return s; 02243 } 02244 02245 //@} 02246 02247 } // namespace std 02248 02249 namespace vigra { 02250 02251 namespace linalg { 02252 02253 namespace detail { 02254 02255 template <class T1, class C1, class T2, class C2, class T3, class C3> 02256 void 02257 columnStatisticsImpl(MultiArrayView<2, T1, C1> const & A, 02258 MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & sumOfSquaredDifferences) 02259 { 02260 MultiArrayIndex m = rowCount(A); 02261 MultiArrayIndex n = columnCount(A); 02262 vigra_precondition(1 == rowCount(mean) && n == columnCount(mean) && 02263 1 == rowCount(sumOfSquaredDifferences) && n == columnCount(sumOfSquaredDifferences), 02264 "columnStatistics(): Shape mismatch between input and output."); 02265 02266 // West's algorithm for incremental variance computation 02267 mean.init(NumericTraits<T2>::zero()); 02268 sumOfSquaredDifferences.init(NumericTraits<T3>::zero()); 02269 02270 for(MultiArrayIndex k=0; k<m; ++k) 02271 { 02272 Matrix<T2> t = rowVector(A, k) - mean; 02273 typename NumericTraits<T2>::RealPromote f = 1.0 / (k + 1.0), 02274 f1 = 1.0 - f; 02275 mean += f*t; 02276 sumOfSquaredDifferences += f1*sq(t); 02277 } 02278 } 02279 02280 template <class T1, class C1, class T2, class C2, class T3, class C3> 02281 void 02282 columnStatistics2PassImpl(MultiArrayView<2, T1, C1> const & A, 02283 MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & sumOfSquaredDifferences) 02284 { 02285 MultiArrayIndex m = rowCount(A); 02286 MultiArrayIndex n = columnCount(A); 02287 vigra_precondition(1 == rowCount(mean) && n == columnCount(mean) && 02288 1 == rowCount(sumOfSquaredDifferences) && n == columnCount(sumOfSquaredDifferences), 02289 "columnStatistics(): Shape mismatch between input and output."); 02290 02291 // two-pass algorithm for incremental variance computation 02292 mean.init(NumericTraits<T2>::zero()); 02293 for(MultiArrayIndex k=0; k<m; ++k) 02294 { 02295 mean += rowVector(A, k); 02296 } 02297 mean /= (double)m; 02298 02299 sumOfSquaredDifferences.init(NumericTraits<T3>::zero()); 02300 for(MultiArrayIndex k=0; k<m; ++k) 02301 { 02302 sumOfSquaredDifferences += sq(rowVector(A, k) - mean); 02303 } 02304 } 02305 02306 } // namespace detail 02307 02308 /** \addtogroup LinearAlgebraFunctions 02309 */ 02310 //@{ 02311 /** Compute statistics of every column of matrix \a A. 02312 02313 The result matrices must be row vectors with as many columns as \a A. 02314 02315 <b> Declarations:</b> 02316 02317 compute only the mean: 02318 \code 02319 namespace vigra { namespace linalg { 02320 template <class T1, class C1, class T2, class C2> 02321 void 02322 columnStatistics(MultiArrayView<2, T1, C1> const & A, 02323 MultiArrayView<2, T2, C2> & mean); 02324 } } 02325 \endcode 02326 02327 compute mean and standard deviation: 02328 \code 02329 namespace vigra { namespace linalg { 02330 template <class T1, class C1, class T2, class C2, class T3, class C3> 02331 void 02332 columnStatistics(MultiArrayView<2, T1, C1> const & A, 02333 MultiArrayView<2, T2, C2> & mean, 02334 MultiArrayView<2, T3, C3> & stdDev); 02335 } } 02336 \endcode 02337 02338 compute mean, standard deviation, and norm: 02339 \code 02340 namespace vigra { namespace linalg { 02341 template <class T1, class C1, class T2, class C2, class T3, class C3, class T4, class C4> 02342 void 02343 columnStatistics(MultiArrayView<2, T1, C1> const & A, 02344 MultiArrayView<2, T2, C2> & mean, 02345 MultiArrayView<2, T3, C3> & stdDev, 02346 MultiArrayView<2, T4, C4> & norm); 02347 } } 02348 \endcode 02349 02350 <b> Usage:</b> 02351 02352 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 02353 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 02354 Namespaces: vigra and vigra::linalg 02355 02356 \code 02357 Matrix A(rows, columns); 02358 .. // fill A 02359 Matrix columnMean(1, columns), columnStdDev(1, columns), columnNorm(1, columns); 02360 02361 columnStatistics(A, columnMean, columnStdDev, columnNorm); 02362 02363 \endcode 02364 */ 02365 doxygen_overloaded_function(template <...> void columnStatistics) 02366 02367 template <class T1, class C1, class T2, class C2> 02368 void 02369 columnStatistics(MultiArrayView<2, T1, C1> const & A, 02370 MultiArrayView<2, T2, C2> & mean) 02371 { 02372 MultiArrayIndex m = rowCount(A); 02373 MultiArrayIndex n = columnCount(A); 02374 vigra_precondition(1 == rowCount(mean) && n == columnCount(mean), 02375 "columnStatistics(): Shape mismatch between input and output."); 02376 02377 mean.init(NumericTraits<T2>::zero()); 02378 02379 for(MultiArrayIndex k=0; k<m; ++k) 02380 { 02381 mean += rowVector(A, k); 02382 } 02383 mean /= T2(m); 02384 } 02385 02386 template <class T1, class C1, class T2, class C2, class T3, class C3> 02387 void 02388 columnStatistics(MultiArrayView<2, T1, C1> const & A, 02389 MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & stdDev) 02390 { 02391 detail::columnStatisticsImpl(A, mean, stdDev); 02392 02393 if(rowCount(A) > 1) 02394 stdDev = sqrt(stdDev / T3(rowCount(A) - 1.0)); 02395 } 02396 02397 template <class T1, class C1, class T2, class C2, class T3, class C3, class T4, class C4> 02398 void 02399 columnStatistics(MultiArrayView<2, T1, C1> const & A, 02400 MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & stdDev, MultiArrayView<2, T4, C4> & norm) 02401 { 02402 MultiArrayIndex m = rowCount(A); 02403 MultiArrayIndex n = columnCount(A); 02404 vigra_precondition(1 == rowCount(mean) && n == columnCount(mean) && 02405 1 == rowCount(stdDev) && n == columnCount(stdDev) && 02406 1 == rowCount(norm) && n == columnCount(norm), 02407 "columnStatistics(): Shape mismatch between input and output."); 02408 02409 detail::columnStatisticsImpl(A, mean, stdDev); 02410 norm = sqrt(stdDev + T2(m) * sq(mean)); 02411 stdDev = sqrt(stdDev / T3(m - 1.0)); 02412 } 02413 02414 /** Compute statistics of every row of matrix \a A. 02415 02416 The result matrices must be column vectors with as many rows as \a A. 02417 02418 <b> Declarations:</b> 02419 02420 compute only the mean: 02421 \code 02422 namespace vigra { namespace linalg { 02423 template <class T1, class C1, class T2, class C2> 02424 void 02425 rowStatistics(MultiArrayView<2, T1, C1> const & A, 02426 MultiArrayView<2, T2, C2> & mean); 02427 } } 02428 \endcode 02429 02430 compute mean and standard deviation: 02431 \code 02432 namespace vigra { namespace linalg { 02433 template <class T1, class C1, class T2, class C2, class T3, class C3> 02434 void 02435 rowStatistics(MultiArrayView<2, T1, C1> const & A, 02436 MultiArrayView<2, T2, C2> & mean, 02437 MultiArrayView<2, T3, C3> & stdDev); 02438 } } 02439 \endcode 02440 02441 compute mean, standard deviation, and norm: 02442 \code 02443 namespace vigra { namespace linalg { 02444 template <class T1, class C1, class T2, class C2, class T3, class C3, class T4, class C4> 02445 void 02446 rowStatistics(MultiArrayView<2, T1, C1> const & A, 02447 MultiArrayView<2, T2, C2> & mean, 02448 MultiArrayView<2, T3, C3> & stdDev, 02449 MultiArrayView<2, T4, C4> & norm); 02450 } } 02451 \endcode 02452 02453 <b> Usage:</b> 02454 02455 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 02456 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 02457 Namespaces: vigra and vigra::linalg 02458 02459 \code 02460 Matrix A(rows, columns); 02461 .. // fill A 02462 Matrix rowMean(rows, 1), rowStdDev(rows, 1), rowNorm(rows, 1); 02463 02464 rowStatistics(a, rowMean, rowStdDev, rowNorm); 02465 02466 \endcode 02467 */ 02468 doxygen_overloaded_function(template <...> void rowStatistics) 02469 02470 template <class T1, class C1, class T2, class C2> 02471 void 02472 rowStatistics(MultiArrayView<2, T1, C1> const & A, 02473 MultiArrayView<2, T2, C2> & mean) 02474 { 02475 vigra_precondition(1 == columnCount(mean) && rowCount(A) == rowCount(mean), 02476 "rowStatistics(): Shape mismatch between input and output."); 02477 MultiArrayView<2, T2, StridedArrayTag> tm = transpose(mean); 02478 columnStatistics(transpose(A), tm); 02479 } 02480 02481 template <class T1, class C1, class T2, class C2, class T3, class C3> 02482 void 02483 rowStatistics(MultiArrayView<2, T1, C1> const & A, 02484 MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & stdDev) 02485 { 02486 vigra_precondition(1 == columnCount(mean) && rowCount(A) == rowCount(mean) && 02487 1 == columnCount(stdDev) && rowCount(A) == rowCount(stdDev), 02488 "rowStatistics(): Shape mismatch between input and output."); 02489 MultiArrayView<2, T2, StridedArrayTag> tm = transpose(mean); 02490 MultiArrayView<2, T3, StridedArrayTag> ts = transpose(stdDev); 02491 columnStatistics(transpose(A), tm, ts); 02492 } 02493 02494 template <class T1, class C1, class T2, class C2, class T3, class C3, class T4, class C4> 02495 void 02496 rowStatistics(MultiArrayView<2, T1, C1> const & A, 02497 MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & stdDev, MultiArrayView<2, T4, C4> & norm) 02498 { 02499 vigra_precondition(1 == columnCount(mean) && rowCount(A) == rowCount(mean) && 02500 1 == columnCount(stdDev) && rowCount(A) == rowCount(stdDev) && 02501 1 == columnCount(norm) && rowCount(A) == rowCount(norm), 02502 "rowStatistics(): Shape mismatch between input and output."); 02503 MultiArrayView<2, T2, StridedArrayTag> tm = transpose(mean); 02504 MultiArrayView<2, T3, StridedArrayTag> ts = transpose(stdDev); 02505 MultiArrayView<2, T4, StridedArrayTag> tn = transpose(norm); 02506 columnStatistics(transpose(A), tm, ts, tn); 02507 } 02508 02509 namespace detail { 02510 02511 template <class T1, class C1, class U, class T2, class C2, class T3, class C3> 02512 void updateCovarianceMatrix(MultiArrayView<2, T1, C1> const & features, 02513 U & count, MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & covariance) 02514 { 02515 MultiArrayIndex n = std::max(rowCount(features), columnCount(features)); 02516 vigra_precondition(std::min(rowCount(features), columnCount(features)) == 1, 02517 "updateCovarianceMatrix(): Features must be a row or column vector."); 02518 vigra_precondition(mean.shape() == features.shape(), 02519 "updateCovarianceMatrix(): Shape mismatch between feature vector and mean vector."); 02520 vigra_precondition(n == rowCount(covariance) && n == columnCount(covariance), 02521 "updateCovarianceMatrix(): Shape mismatch between feature vector and covariance matrix."); 02522 02523 // West's algorithm for incremental covariance matrix computation 02524 Matrix<T2> t = features - mean; 02525 ++count; 02526 double f = 1.0 / count, 02527 f1 = 1.0 - f; 02528 mean += f*t; 02529 02530 if(rowCount(features) == 1) // update column covariance from current row 02531 { 02532 for(MultiArrayIndex k=0; k<n; ++k) 02533 { 02534 covariance(k, k) += f1*sq(t(0, k)); 02535 for(MultiArrayIndex l=k+1; l<n; ++l) 02536 { 02537 covariance(k, l) += f1*t(0, k)*t(0, l); 02538 covariance(l, k) = covariance(k, l); 02539 } 02540 } 02541 } 02542 else // update row covariance from current column 02543 { 02544 for(MultiArrayIndex k=0; k<n; ++k) 02545 { 02546 covariance(k, k) += f1*sq(t(k, 0)); 02547 for(MultiArrayIndex l=k+1; l<n; ++l) 02548 { 02549 covariance(k, l) += f1*t(k, 0)*t(l, 0); 02550 covariance(l, k) = covariance(k, l); 02551 } 02552 } 02553 } 02554 } 02555 02556 } // namespace detail 02557 02558 /*! Compute the covariance matrix between the columns of a matrix \a features. 02559 02560 The result matrix \a covariance must by a square matrix with as many rows and 02561 columns as the number of columns in matrix \a features. 02562 02563 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>><br> 02564 Namespace: vigra 02565 */ 02566 template <class T1, class C1, class T2, class C2> 02567 void covarianceMatrixOfColumns(MultiArrayView<2, T1, C1> const & features, 02568 MultiArrayView<2, T2, C2> & covariance) 02569 { 02570 MultiArrayIndex m = rowCount(features), n = columnCount(features); 02571 vigra_precondition(n == rowCount(covariance) && n == columnCount(covariance), 02572 "covarianceMatrixOfColumns(): Shape mismatch between feature matrix and covariance matrix."); 02573 MultiArrayIndex count = 0; 02574 Matrix<T2> means(1, n); 02575 covariance.init(NumericTraits<T2>::zero()); 02576 for(MultiArrayIndex k=0; k<m; ++k) 02577 detail::updateCovarianceMatrix(rowVector(features, k), count, means, covariance); 02578 covariance /= T2(m - 1); 02579 } 02580 02581 /*! Compute the covariance matrix between the columns of a matrix \a features. 02582 02583 The result is returned as a square temporary matrix with as many rows and 02584 columns as the number of columns in matrix \a features. 02585 02586 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>><br> 02587 Namespace: vigra 02588 */ 02589 template <class T, class C> 02590 TemporaryMatrix<T> 02591 covarianceMatrixOfColumns(MultiArrayView<2, T, C> const & features) 02592 { 02593 TemporaryMatrix<T> res(columnCount(features), columnCount(features)); 02594 covarianceMatrixOfColumns(features, res); 02595 return res; 02596 } 02597 02598 /*! Compute the covariance matrix between the rows of a matrix \a features. 02599 02600 The result matrix \a covariance must by a square matrix with as many rows and 02601 columns as the number of rows in matrix \a features. 02602 02603 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>><br> 02604 Namespace: vigra 02605 */ 02606 template <class T1, class C1, class T2, class C2> 02607 void covarianceMatrixOfRows(MultiArrayView<2, T1, C1> const & features, 02608 MultiArrayView<2, T2, C2> & covariance) 02609 { 02610 MultiArrayIndex m = rowCount(features), n = columnCount(features); 02611 vigra_precondition(m == rowCount(covariance) && m == columnCount(covariance), 02612 "covarianceMatrixOfRows(): Shape mismatch between feature matrix and covariance matrix."); 02613 MultiArrayIndex count = 0; 02614 Matrix<T2> means(m, 1); 02615 covariance.init(NumericTraits<T2>::zero()); 02616 for(MultiArrayIndex k=0; k<n; ++k) 02617 detail::updateCovarianceMatrix(columnVector(features, k), count, means, covariance); 02618 covariance /= T2(m - 1); 02619 } 02620 02621 /*! Compute the covariance matrix between the rows of a matrix \a features. 02622 02623 The result is returned as a square temporary matrix with as many rows and 02624 columns as the number of rows in matrix \a features. 02625 02626 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>><br> 02627 Namespace: vigra 02628 */ 02629 template <class T, class C> 02630 TemporaryMatrix<T> 02631 covarianceMatrixOfRows(MultiArrayView<2, T, C> const & features) 02632 { 02633 TemporaryMatrix<T> res(rowCount(features), rowCount(features)); 02634 covarianceMatrixOfRows(features, res); 02635 return res; 02636 } 02637 02638 enum DataPreparationGoals { ZeroMean = 1, UnitVariance = 2, UnitNorm = 4 }; 02639 02640 inline DataPreparationGoals operator|(DataPreparationGoals l, DataPreparationGoals r) 02641 { 02642 return DataPreparationGoals(int(l) | int(r)); 02643 } 02644 02645 namespace detail { 02646 02647 template <class T, class C1, class C2, class C3, class C4> 02648 void 02649 prepareDataImpl(const MultiArrayView<2, T, C1> & A, 02650 MultiArrayView<2, T, C2> & res, MultiArrayView<2, T, C3> & offset, MultiArrayView<2, T, C4> & scaling, 02651 DataPreparationGoals goals) 02652 { 02653 MultiArrayIndex m = rowCount(A); 02654 MultiArrayIndex n = columnCount(A); 02655 vigra_precondition(A.shape() == res.shape() && 02656 n == columnCount(offset) && 1 == rowCount(offset) && 02657 n == columnCount(scaling) && 1 == rowCount(scaling), 02658 "prepareDataImpl(): Shape mismatch between input and output."); 02659 02660 if(!goals) 02661 { 02662 res = A; 02663 offset.init(NumericTraits<T>::zero()); 02664 scaling.init(NumericTraits<T>::one()); 02665 return; 02666 } 02667 02668 bool zeroMean = (goals & ZeroMean) != 0; 02669 bool unitVariance = (goals & UnitVariance) != 0; 02670 bool unitNorm = (goals & UnitNorm) != 0; 02671 02672 vigra_precondition(!(unitVariance && unitNorm), 02673 "prepareDataImpl(): Unit variance and unit norm cannot be achieved at the same time."); 02674 02675 Matrix<T> mean(1, n), sumOfSquaredDifferences(1, n); 02676 detail::columnStatisticsImpl(A, mean, sumOfSquaredDifferences); 02677 02678 for(MultiArrayIndex k=0; k<n; ++k) 02679 { 02680 T stdDev = std::sqrt(sumOfSquaredDifferences(0, k) / T(m-1)); 02681 if(closeAtTolerance(stdDev / mean(0,k), NumericTraits<T>::zero())) 02682 stdDev = NumericTraits<T>::zero(); 02683 if(zeroMean && stdDev > NumericTraits<T>::zero()) 02684 { 02685 columnVector(res, k) = columnVector(A, k) - mean(0,k); 02686 offset(0, k) = mean(0, k); 02687 mean(0, k) = NumericTraits<T>::zero(); 02688 } 02689 else 02690 { 02691 columnVector(res, k) = columnVector(A, k); 02692 offset(0, k) = NumericTraits<T>::zero(); 02693 } 02694 02695 T norm = mean(0,k) == NumericTraits<T>::zero() 02696 ? std::sqrt(sumOfSquaredDifferences(0, k)) 02697 : std::sqrt(sumOfSquaredDifferences(0, k) + T(m) * sq(mean(0,k))); 02698 if(unitNorm && norm > NumericTraits<T>::zero()) 02699 { 02700 columnVector(res, k) /= norm; 02701 scaling(0, k) = NumericTraits<T>::one() / norm; 02702 } 02703 else if(unitVariance && stdDev > NumericTraits<T>::zero()) 02704 { 02705 columnVector(res, k) /= stdDev; 02706 scaling(0, k) = NumericTraits<T>::one() / stdDev; 02707 } 02708 else 02709 { 02710 scaling(0, k) = NumericTraits<T>::one(); 02711 } 02712 } 02713 } 02714 02715 } // namespace detail 02716 02717 /*! Standardize the columns of a matrix according to given <tt>DataPreparationGoals</tt>. 02718 02719 For every column of the matrix \a A, this function computes mean, 02720 standard deviation, and norm. It then applies a linear transformation to the values of 02721 the column according to these statistics and the given <tt>DataPreparationGoals</tt>. 02722 The result is returned in matrix \a res which must have the same size as \a A. 02723 Optionally, the transformation applied can also be returned in the matrices \a offset 02724 and \a scaling (see below for an example how these matrices can be used to standardize 02725 more data according to the same transformation). 02726 02727 The following <tt>DataPreparationGoals</tt> are supported: 02728 02729 <DL> 02730 <DT><tt>ZeroMean</tt><DD> Subtract the column mean form every column if the values in the column are not constant. 02731 Do nothing in a constant column. 02732 <DT><tt>UnitVariance</tt><DD> Divide by the column standard deviation if the values in the column are not constant. 02733 Do nothing in a constant column. 02734 <DT><tt>UnitNorm</tt><DD> Divide by the column norm if it is non-zero. 02735 <DT><tt>ZeroMean | UnitVariance</tt><DD> First subtact the mean and then divide by the standard deviation, unless the 02736 column is constant (in which case the column remains unchanged). 02737 <DT><tt>ZeroMean | UnitNorm</tt><DD> If the column is non-constant, subtract the mean. Then divide by the norm 02738 of the result if the norm is non-zero. 02739 </DL> 02740 02741 <b> Declarations:</b> 02742 02743 Standardize the matrix and return the parameters of the linear transformation. 02744 The matrices \a offset and \a scaling must be row vectors with as many columns as \a A. 02745 \code 02746 namespace vigra { namespace linalg { 02747 template <class T, class C1, class C2, class C3, class C4> 02748 void 02749 prepareColumns(MultiArrayView<2, T, C1> const & A, 02750 MultiArrayView<2, T, C2> & res, 02751 MultiArrayView<2, T, C3> & offset, 02752 MultiArrayView<2, T, C4> & scaling, 02753 DataPreparationGoals goals = ZeroMean | UnitVariance); 02754 } } 02755 \endcode 02756 02757 Only standardize the matrix. 02758 \code 02759 namespace vigra { namespace linalg { 02760 template <class T, class C1, class C2> 02761 void 02762 prepareColumns(MultiArrayView<2, T, C1> const & A, 02763 MultiArrayView<2, T, C2> & res, 02764 DataPreparationGoals goals = ZeroMean | UnitVariance); 02765 } } 02766 \endcode 02767 02768 <b> Usage:</b> 02769 02770 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 02771 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 02772 Namespaces: vigra and vigra::linalg 02773 02774 \code 02775 Matrix A(rows, columns); 02776 .. // fill A 02777 Matrix standardizedA(rows, columns), offset(1, columns), scaling(1, columns); 02778 02779 prepareColumns(A, standardizedA, offset, scaling, ZeroMean | UnitNorm); 02780 02781 // use offset and scaling to prepare additional data according to the same transformation 02782 Matrix newData(nrows, columns); 02783 02784 Matrix standardizedNewData = (newData - repeatMatrix(offset, nrows, 1)) * pointWise(repeatMatrix(scaling, nrows, 1)); 02785 02786 \endcode 02787 */ 02788 doxygen_overloaded_function(template <...> void prepareColumns) 02789 02790 template <class T, class C1, class C2, class C3, class C4> 02791 inline void 02792 prepareColumns(MultiArrayView<2, T, C1> const & A, 02793 MultiArrayView<2, T, C2> & res, MultiArrayView<2, T, C3> & offset, MultiArrayView<2, T, C4> & scaling, 02794 DataPreparationGoals goals = ZeroMean | UnitVariance) 02795 { 02796 detail::prepareDataImpl(A, res, offset, scaling, goals); 02797 } 02798 02799 template <class T, class C1, class C2> 02800 inline void 02801 prepareColumns(MultiArrayView<2, T, C1> const & A, MultiArrayView<2, T, C2> & res, 02802 DataPreparationGoals goals = ZeroMean | UnitVariance) 02803 { 02804 Matrix<T> offset(1, columnCount(A)), scaling(1, columnCount(A)); 02805 detail::prepareDataImpl(A, res, offset, scaling, goals); 02806 } 02807 02808 /*! Standardize the rows of a matrix according to given <tt>DataPreparationGoals</tt>. 02809 02810 This algorithm works in the same way as \ref prepareColumns() (see there for detailed 02811 documentation), but is applied to the rows of the matrix \a A instead. Accordingly, the 02812 matrices holding the parameters of the linear transformation must be column vectors 02813 with as many rows as \a A. 02814 02815 <b> Declarations:</b> 02816 02817 Standardize the matrix and return the parameters of the linear transformation. 02818 The matrices \a offset and \a scaling must be column vectors 02819 with as many rows as \a A. 02820 02821 \code 02822 namespace vigra { namespace linalg { 02823 template <class T, class C1, class C2, class C3, class C4> 02824 void 02825 prepareRows(MultiArrayView<2, T, C1> const & A, 02826 MultiArrayView<2, T, C2> & res, 02827 MultiArrayView<2, T, C3> & offset, 02828 MultiArrayView<2, T, C4> & scaling, 02829 DataPreparationGoals goals = ZeroMean | UnitVariance)´; 02830 } } 02831 \endcode 02832 02833 Only standardize the matrix. 02834 \code 02835 namespace vigra { namespace linalg { 02836 template <class T, class C1, class C2> 02837 void 02838 prepareRows(MultiArrayView<2, T, C1> const & A, 02839 MultiArrayView<2, T, C2> & res, 02840 DataPreparationGoals goals = ZeroMean | UnitVariance); 02841 } } 02842 \endcode 02843 02844 <b> Usage:</b> 02845 02846 <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br> 02847 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 02848 Namespaces: vigra and vigra::linalg 02849 02850 \code 02851 Matrix A(rows, columns); 02852 .. // fill A 02853 Matrix standardizedA(rows, columns), offset(rows, 1), scaling(rows, 1); 02854 02855 prepareRows(A, standardizedA, offset, scaling, ZeroMean | UnitNorm); 02856 02857 // use offset and scaling to prepare additional data according to the same transformation 02858 Matrix newData(rows, ncolumns); 02859 02860 Matrix standardizedNewData = (newData - repeatMatrix(offset, 1, ncolumns)) * pointWise(repeatMatrix(scaling, 1, ncolumns)); 02861 02862 \endcode 02863 */ 02864 doxygen_overloaded_function(template <...> void prepareRows) 02865 02866 template <class T, class C1, class C2, class C3, class C4> 02867 inline void 02868 prepareRows(MultiArrayView<2, T, C1> const & A, 02869 MultiArrayView<2, T, C2> & res, MultiArrayView<2, T, C3> & offset, MultiArrayView<2, T, C4> & scaling, 02870 DataPreparationGoals goals = ZeroMean | UnitVariance) 02871 { 02872 MultiArrayView<2, T, StridedArrayTag> tr = transpose(res), to = transpose(offset), ts = transpose(scaling); 02873 detail::prepareDataImpl(transpose(A), tr, to, ts, goals); 02874 } 02875 02876 template <class T, class C1, class C2> 02877 inline void 02878 prepareRows(MultiArrayView<2, T, C1> const & A, MultiArrayView<2, T, C2> & res, 02879 DataPreparationGoals goals = ZeroMean | UnitVariance) 02880 { 02881 MultiArrayView<2, T, StridedArrayTag> tr = transpose(res); 02882 Matrix<T> offset(rowCount(A), 1), scaling(rowCount(A), 1); 02883 detail::prepareDataImpl(transpose(A), tr, offset, scaling, goals); 02884 } 02885 02886 //@} 02887 02888 } // namespace linalg 02889 02890 using linalg::columnStatistics; 02891 using linalg::prepareColumns; 02892 using linalg::rowStatistics; 02893 using linalg::prepareRows; 02894 using linalg::ZeroMean; 02895 using linalg::UnitVariance; 02896 using linalg::UnitNorm; 02897 02898 } // namespace vigra 02899 02900 02901 02902 #endif // VIGRA_MATRIX_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|