[ 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)