[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

vigra/matrix.hxx |

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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.7.1 (3 Dec 2010)