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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.7.0 (Thu Aug 25 2011)