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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python VIGRA 1.6.0 (13 Aug 2008)