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