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

vigra/matrix.hxx |

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