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

vigra/matrix.hxx VIGRA

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