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