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

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

html generated using doxygen and Python
VIGRA 1.6.0 (5 Nov 2009)