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