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