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

vigra/matrix.hxx VIGRA

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*     Copyright 2003-2008 by Gunnar Kedenburg and Ullrich Koethe       */
00004 /*                                                                      */
00005 /*    This file is part of the VIGRA computer vision library.           */
00006 /*    The VIGRA Website is                                              */
00007 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
00008 /*    Please direct questions, bug reports, and contributions to        */
00009 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00010 /*        vigra@informatik.uni-hamburg.de                               */
00011 /*                                                                      */
00012 /*    Permission is hereby granted, free of charge, to any person       */
00013 /*    obtaining a copy of this software and associated documentation    */
00014 /*    files (the "Software"), to deal in the Software without           */
00015 /*    restriction, including without limitation the rights to use,      */
00016 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00017 /*    sell copies of the Software, and to permit persons to whom the    */
00018 /*    Software is furnished to do so, subject to the following          */
00019 /*    conditions:                                                       */
00020 /*                                                                      */
00021 /*    The above copyright notice and this permission notice shall be    */
00022 /*    included in all copies or substantial portions of the             */
00023 /*    Software.                                                         */
00024 /*                                                                      */
00025 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00026 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00027 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00028 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00029 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00030 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00031 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00032 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
00033 /*                                                                      */
00034 /************************************************************************/
00035 
00036 
00037 #ifndef VIGRA_MATRIX_HXX
00038 #define VIGRA_MATRIX_HXX
00039 
00040 #include <cmath>
00041 #include <iosfwd>
00042 #include <iomanip>
00043 #include "multi_array.hxx"
00044 #include "mathutil.hxx"
00045 #include "numerictraits.hxx"
00046 #include "multi_pointoperators.hxx"
00047 
00048 
00049 namespace vigra
00050 {
00051 
00052 /** \defgroup LinearAlgebraModule Linear Algebra
00053 
00054     \brief Classes and functions for matrix algebra, linear equations systems, eigen systems, least squares etc.
00055 */
00056 
00057 /** \ingroup LinearAlgebraModule
00058 
00059     Namespace <tt>vigra/linalg</tt> hold VIGRA's linear algebra functionality. But most of its contents
00060     is exported into namespace <tt>vigra</tt> via <tt>using</tt> directives.
00061 */
00062 namespace linalg
00063 {
00064 
00065 template <class T, class C>
00066 inline MultiArrayIndex 
00067 rowCount(const MultiArrayView<2, T, C> &x);
00068 
00069 template <class T, class C>
00070 inline MultiArrayIndex 
00071 columnCount(const MultiArrayView<2, T, C> &x);
00072 
00073 template <class T, class C>
00074 inline MultiArrayView <2, T, C>
00075 rowVector(MultiArrayView <2, T, C> const & m, MultiArrayIndex d);
00076 
00077 template <class T, class C>
00078 inline MultiArrayView <2, T, C>
00079 columnVector(MultiArrayView<2, T, C> const & m, MultiArrayIndex d);
00080 
00081 template <class T, class ALLOC = std::allocator<T> >
00082 class TemporaryMatrix;
00083 
00084 template <class T, class C1, class C2>
00085 void transpose(const MultiArrayView<2, T, C1> &v, MultiArrayView<2, T, C2> &r);
00086 
00087 template <class T, class C>
00088 bool isSymmetric(const MultiArrayView<2, T, C> &v);
00089 
00090 enum RawArrayMemoryLayout { RowMajor, ColumnMajor };
00091 
00092 /********************************************************/
00093 /*                                                      */
00094 /*                        Matrix                        */
00095 /*                                                      */
00096 /********************************************************/
00097 
00098 /** Matrix class. 
00099 
00100     \ingroup LinearAlgebraModule 
00101 
00102     This is the basic class for all linear algebra computations. Matrices are
00103     stored in a <i>column-major</i> format, i.e. the row index is varying fastest.
00104     This is the same format as in the lapack and gmm++ libraries, so it will
00105     be easy to interface these libraries. In fact, if you need optimized
00106     high performance code, you should use them. The VIGRA linear algebra
00107     functionality is provided for smaller problems and rapid prototyping
00108     (no one wants to spend half the day installing a new library just to
00109     discover that the new algorithm idea didn't work anyway).
00110 
00111     <b>See also:</b>
00112     <ul>
00113     <li> \ref LinearAlgebraFunctions
00114     </ul>
00115 
00116     <b>\#include</b> <vigra/matrix.hxx> or<br>
00117     <b>\#include</b> <vigra/linear_algebra.hxx><br>
00118         Namespaces: vigra and vigra::linalg
00119 */
00120 template <class T, class ALLOC = std::allocator<T> >
00121 class Matrix
00122 : public MultiArray<2, T, ALLOC>
00123 {
00124     typedef MultiArray<2, T, ALLOC> BaseType;
00125 
00126   public:
00127     typedef Matrix<T, ALLOC>                        matrix_type;
00128     typedef TemporaryMatrix<T, ALLOC>               temp_type;
00129     typedef MultiArrayView<2, T, UnstridedArrayTag> view_type;
00130     typedef typename BaseType::value_type           value_type;
00131     typedef typename BaseType::pointer              pointer;
00132     typedef typename BaseType::const_pointer        const_pointer;
00133     typedef typename BaseType::reference            reference;
00134     typedef typename BaseType::const_reference      const_reference;
00135     typedef typename BaseType::difference_type      difference_type;
00136     typedef typename BaseType::difference_type_1    difference_type_1;
00137     typedef ALLOC                                   allocator_type;
00138 
00139         /** default constructor
00140          */
00141     Matrix()
00142     {}
00143 
00144         /** construct with given allocator
00145          */
00146     explicit Matrix(ALLOC const & alloc)
00147     : BaseType(alloc)
00148     {}
00149 
00150         /** construct with given shape and init all
00151             elements with zero. Note that the order of the axes is
00152             <tt>difference_type(rows, columns)</tt> which
00153             is the opposite of the usual VIGRA convention.
00154          */
00155     explicit Matrix(const difference_type &shape,
00156                     ALLOC const & alloc = allocator_type())
00157     : BaseType(shape, alloc)
00158     {}
00159 
00160         /** construct with given shape and init all
00161             elements with zero. Note that the order of the axes is
00162             <tt>(rows, columns)</tt> which
00163             is the opposite of the usual VIGRA convention.
00164          */
00165     Matrix(difference_type_1 rows, difference_type_1 columns,
00166                     ALLOC const & alloc = allocator_type())
00167     : BaseType(difference_type(rows, columns), alloc)
00168     {}
00169 
00170         /** construct with given shape and init all
00171             elements with the constant \a init. Note that the order of the axes is
00172             <tt>difference_type(rows, columns)</tt> which
00173             is the opposite of the usual VIGRA convention.
00174          */
00175     Matrix(const difference_type &shape, const_reference init,
00176            allocator_type const & alloc = allocator_type())
00177     : BaseType(shape, init, alloc)
00178     {}
00179 
00180         /** construct with given shape and init all
00181             elements with the constant \a init. Note that the order of the axes is
00182             <tt>(rows, columns)</tt> which
00183             is the opposite of the usual VIGRA convention.
00184          */
00185     Matrix(difference_type_1 rows, difference_type_1 columns, const_reference init,
00186            allocator_type const & alloc = allocator_type())
00187     : BaseType(difference_type(rows, columns), init, alloc)
00188     {}
00189 
00190         /** construct with given shape and copy data from C-style array \a init.
00191             Unless \a layout is <tt>ColumnMajor</tt>, the elements in this array
00192             are assumed to be given in row-major order (the C standard order) and
00193             will automatically be converted to the required column-major format.
00194             Note that the order of the axes is <tt>difference_type(rows, columns)</tt> which
00195             is the opposite of the usual VIGRA convention.
00196          */
00197     Matrix(const difference_type &shape, const_pointer init, RawArrayMemoryLayout layout = RowMajor,
00198            allocator_type const & alloc = allocator_type())
00199     : BaseType(shape, alloc) // FIXME: this function initializes the memory twice
00200     {
00201         if(layout == RowMajor)
00202         {
00203             difference_type trans(shape[1], shape[0]);
00204             linalg::transpose(MultiArrayView<2, T>(trans, const_cast<pointer>(init)), *this);
00205         }
00206         else
00207         {
00208             std::copy(init, init + elementCount(), this->data());
00209         }
00210     }
00211 
00212         /** construct with given shape and copy data from C-style array \a init.
00213             Unless \a layout is <tt>ColumnMajor</tt>, the elements in this array
00214             are assumed to be given in row-major order (the C standard order) and
00215             will automatically be converted to the required column-major format.
00216             Note that the order of the axes is <tt>(rows, columns)</tt> which
00217             is the opposite of the usual VIGRA convention.
00218          */
00219     Matrix(difference_type_1 rows, difference_type_1 columns, const_pointer init, RawArrayMemoryLayout layout = RowMajor,
00220            allocator_type const & alloc = allocator_type())
00221     : BaseType(difference_type(rows, columns), alloc) // FIXME: this function initializes the memory twice
00222     {
00223         if(layout == RowMajor)
00224         {
00225             difference_type trans(columns, rows);
00226             linalg::transpose(MultiArrayView<2, T>(trans, const_cast<pointer>(init)), *this);
00227         }
00228         else
00229         {
00230             std::copy(init, init + elementCount(), this->data());
00231         }
00232     }
00233 
00234         /** copy constructor. Allocates new memory and
00235             copies tha data.
00236          */
00237     Matrix(const Matrix &rhs)
00238     : BaseType(rhs)
00239     {}
00240 
00241         /** construct from temporary matrix, which looses its data.
00242 
00243             This operation is equivalent to
00244             \code
00245                 TemporaryMatrix<T> temp = ...;
00246 
00247                 Matrix<T> m;
00248                 m.swap(temp);
00249             \endcode
00250          */
00251     Matrix(const TemporaryMatrix<T, ALLOC> &rhs)
00252     : BaseType(rhs.allocator())
00253     {
00254         this->swap(const_cast<TemporaryMatrix<T, ALLOC> &>(rhs));
00255     }
00256 
00257         /** construct from a MultiArrayView. Allocates new memory and
00258             copies tha data. \a rhs is assumed to be in column-major order already.
00259          */
00260     template<class U, class C>
00261     Matrix(const MultiArrayView<2, U, C> &rhs)
00262     : BaseType(rhs)
00263     {}
00264 
00265         /** assignment.
00266             If the size of \a rhs is the same as the matrix's old size, only the data
00267             are copied. Otherwise, new storage is allocated, which invalidates
00268             all objects (array views, iterators) depending on the matrix.
00269          */
00270     Matrix & operator=(const Matrix &rhs)
00271     {
00272         BaseType::operator=(rhs); // has the correct semantics already
00273         return *this;
00274     }
00275 
00276         /** assign a temporary matrix. If the shapes of the two matrices match,
00277             only the data are copied (in order to not invalidate views and iterators
00278             depending on this matrix). Otherwise, the memory is swapped
00279             between the two matrices, so that all depending objects
00280             (array views, iterators) ar invalidated.
00281          */
00282     Matrix & operator=(const TemporaryMatrix<T, ALLOC> &rhs)
00283     {
00284         if(this->shape() == rhs.shape())
00285             this->copy(rhs);
00286         else
00287             this->swap(const_cast<TemporaryMatrix<T, ALLOC> &>(rhs));
00288         return *this;
00289     }
00290 
00291         /** assignment from arbitrary 2-dimensional MultiArrayView.<br>
00292             If the size of \a rhs is the same as the matrix's old size, only the data
00293             are copied. Otherwise, new storage is allocated, which invalidates
00294             all objects (array views, iterators) depending on the matrix.
00295             \a rhs is assumed to be in column-major order already.
00296          */
00297     template <class U, class C>
00298     Matrix & operator=(const MultiArrayView<2, U, C> &rhs)
00299     {
00300         BaseType::operator=(rhs); // has the correct semantics already
00301         return *this;
00302     }
00303 
00304         /** 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)