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

matrix.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 2003-2008 by Gunnar Kedenburg and Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 
37 #ifndef VIGRA_MATRIX_HXX
38 #define VIGRA_MATRIX_HXX
39 
40 #include <cmath>
41 #include <iosfwd>
42 #include <iomanip>
43 #include "multi_array.hxx"
44 #include "mathutil.hxx"
45 #include "numerictraits.hxx"
46 
47 
48 namespace vigra
49 {
50 
51 /** \defgroup LinearAlgebraModule Linear Algebra
52 
53  \brief Classes and functions for matrix algebra, linear equations systems, eigen systems, least squares etc.
54 */
55 
56 /** \ingroup LinearAlgebraModule
57 
58  Namespace <tt>vigra/linalg</tt> hold VIGRA's linear algebra functionality. But most of its contents
59  is exported into namespace <tt>vigra</tt> via <tt>using</tt> directives.
60 */
61 namespace linalg
62 {
63 
64 template <class T, class C>
65 inline MultiArrayIndex
66 rowCount(const MultiArrayView<2, T, C> &x);
67 
68 template <class T, class C>
69 inline MultiArrayIndex
70 columnCount(const MultiArrayView<2, T, C> &x);
71 
72 template <class T, class C>
73 inline MultiArrayView <2, T, C>
74 rowVector(MultiArrayView <2, T, C> const & m, MultiArrayIndex d);
75 
76 template <class T, class C>
77 inline MultiArrayView <2, T, C>
78 columnVector(MultiArrayView<2, T, C> const & m, MultiArrayIndex d);
79 
80 template <class T, class ALLOC>
81 class TemporaryMatrix;
82 
83 template <class T, class C1, class C2>
84 void transpose(const MultiArrayView<2, T, C1> &v, MultiArrayView<2, T, C2> &r);
85 
86 template <class T, class C>
87 bool isSymmetric(const MultiArrayView<2, T, C> &v);
88 
89 enum RawArrayMemoryLayout { RowMajor, ColumnMajor };
90 
91 /********************************************************/
92 /* */
93 /* Matrix */
94 /* */
95 /********************************************************/
96 
97 /** Matrix class.
98 
99  \ingroup LinearAlgebraModule
100 
101  This is the basic class for all linear algebra computations. Matrices are
102  strored in a <i>column-major</i> format, i.e. the row index is varying fastest.
103  This is the same format as in the lapack and gmm++ libraries, so it will
104  be easy to interface these libraries. In fact, if you need optimized
105  high performance code, you should use them. The VIGRA linear algebra
106  functionality is provided for smaller problems and rapid prototyping
107  (no one wants to spend half the day installing a new library just to
108  discover that the new algorithm idea didn't work anyway).
109 
110  <b>See also:</b>
111  <ul>
112  <li> \ref LinearAlgebraFunctions
113  </ul>
114 
115  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
116  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
117  Namespaces: vigra and vigra::linalg
118 */
119 template <class T, class ALLOC = std::allocator<T> >
120 class Matrix
121 : public MultiArray<2, T, ALLOC>
122 {
124 
125  public:
127  typedef TemporaryMatrix<T, ALLOC> temp_type;
130  typedef typename BaseType::pointer pointer;
132  typedef typename BaseType::reference reference;
136  typedef ALLOC allocator_type;
137 
138  /** default constructor
139  */
141  {}
142 
143  /** construct with given allocator
144  */
145  explicit Matrix(ALLOC const & alloc)
146  : BaseType(alloc)
147  {}
148 
149  /** construct with given shape and init all
150  elements with zero. Note that the order of the axes is
151  <tt>difference_type(rows, columns)</tt> which
152  is the opposite of the usual VIGRA convention.
153  */
154  explicit Matrix(const difference_type &shape,
155  ALLOC const & alloc = allocator_type())
156  : BaseType(shape, alloc)
157  {}
158 
159  /** construct with given shape and init all
160  elements with zero. Note that the order of the axes is
161  <tt>(rows, columns)</tt> which
162  is the opposite of the usual VIGRA convention.
163  */
165  ALLOC const & alloc = allocator_type())
166  : BaseType(difference_type(rows, columns), alloc)
167  {}
168 
169  /** construct with given shape and init all
170  elements with the constant \a init. Note that the order of the axes is
171  <tt>difference_type(rows, columns)</tt> which
172  is the opposite of the usual VIGRA convention.
173  */
175  allocator_type const & alloc = allocator_type())
176  : BaseType(shape, init, alloc)
177  {}
178 
179  /** construct with given shape and init all
180  elements with the constant \a init. Note that the order of the axes is
181  <tt>(rows, columns)</tt> which
182  is the opposite of the usual VIGRA convention.
183  */
185  allocator_type const & alloc = allocator_type())
186  : BaseType(difference_type(rows, columns), init, alloc)
187  {}
188 
189  /** construct with given shape and copy data from C-style array \a init.
190  Unless \a layout is <tt>ColumnMajor</tt>, the elements in this array
191  are assumed to be given in row-major order (the C standard order) and
192  will automatically be converted to the required column-major format.
193  Note that the order of the axes is <tt>difference_type(rows, columns)</tt> which
194  is the opposite of the usual VIGRA convention.
195  */
196  Matrix(const difference_type &shape, const_pointer init, RawArrayMemoryLayout layout = RowMajor,
197  allocator_type const & alloc = allocator_type())
198  : BaseType(shape, alloc) // FIXME: this function initializes the memory twice
199  {
200  if(layout == RowMajor)
201  {
202  difference_type trans(shape[1], shape[0]);
203  linalg::transpose(MultiArrayView<2, T>(trans, const_cast<pointer>(init)), *this);
204  }
205  else
206  {
207  std::copy(init, init + elementCount(), this->data());
208  }
209  }
210 
211  /** construct with given shape and copy data from C-style array \a init.
212  Unless \a layout is <tt>ColumnMajor</tt>, the elements in this array
213  are assumed to be given in row-major order (the C standard order) and
214  will automatically be converted to the required column-major format.
215  Note that the order of the axes is <tt>(rows, columns)</tt> which
216  is the opposite of the usual VIGRA convention.
217  */
218  Matrix(difference_type_1 rows, difference_type_1 columns, const_pointer init, RawArrayMemoryLayout layout = RowMajor,
219  allocator_type const & alloc = allocator_type())
220  : BaseType(difference_type(rows, columns), alloc) // FIXME: this function initializes the memory twice
221  {
222  if(layout == RowMajor)
223  {
224  difference_type trans(columns, rows);
225  linalg::transpose(MultiArrayView<2, T>(trans, const_cast<pointer>(init)), *this);
226  }
227  else
228  {
229  std::copy(init, init + elementCount(), this->data());
230  }
231  }
232 
233  /** copy constructor. Allocates new memory and
234  copies tha data.
235  */
236  Matrix(const Matrix &rhs)
237  : BaseType(rhs)
238  {}
239 
240  /** construct from temporary matrix, which looses its data.
241 
242  This operation is equivalent to
243  \code
244  TemporaryMatrix<T> temp = ...;
245 
246  Matrix<T> m;
247  m.swap(temp);
248  \endcode
249  */
250  Matrix(const TemporaryMatrix<T, ALLOC> &rhs)
251  : BaseType(rhs.allocator())
252  {
253  this->swap(const_cast<TemporaryMatrix<T, ALLOC> &>(rhs));
254  }
255 
256  /** construct from a MultiArrayView. Allocates new memory and
257  copies tha data. \a rhs is assumed to be in column-major order already.
258  */
259  template<class U, class C>
261  : BaseType(rhs)
262  {}
263 
264  /** assignment.
265  If the size of \a rhs is the same as the matrix's old size, only the data
266  are copied. Otherwise, new storage is allocated, which invalidates
267  all objects (array views, iterators) depending on the matrix.
268  */
269  Matrix & operator=(const Matrix &rhs)
270  {
271  BaseType::operator=(rhs); // has the correct semantics already
272  return *this;
273  }
274 
275  /** assign a temporary matrix. If the shapes of the two matrices match,
276  only the data are copied (in order to not invalidate views and iterators
277  depending on this matrix). Otherwise, the memory is swapped
278  between the two matrices, so that all depending objects
279  (array views, iterators) ar invalidated.
280  */
281  Matrix & operator=(const TemporaryMatrix<T, ALLOC> &rhs)
282  {
283  if(this->shape() == rhs.shape())
284  this->copy(rhs);
285  else
286  this->swap(const_cast<TemporaryMatrix<T, ALLOC> &>(rhs));
287  return *this;
288  }
289 
290  /** assignment from arbitrary 2-dimensional MultiArrayView.<br>
291  If the size of \a rhs is the same as the matrix's old size, only the data
292  are copied. Otherwise, new storage is allocated, which invalidates
293  all objects (array views, iterators) depending on the matrix.
294  \a rhs is assumed to be in column-major order already.
295  */
296  template <class U, class C>
298  {
299  BaseType::operator=(rhs); // has the correct semantics already
300  return *this;
301  }
302 
303  /** init elements with a constant
304  */
305  template <class U>
306  Matrix & init(const U & init)
307  {
308  BaseType::init(init);
309  return *this;
310  }
311 
312  /** reshape to the given shape and initialize with zero.
313  */
315  {
316  BaseType::reshape(difference_type(rows, columns));
317  }
318 
319  /** reshape to the given shape and initialize with \a init.
320  */
322  {
323  BaseType::reshape(difference_type(rows, columns), init);
324  }
325 
326  /** reshape to the given shape and initialize with zero.
327  */
328  void reshape(difference_type const & shape)
329  {
330  BaseType::reshape(shape);
331  }
332 
333  /** reshape to the given shape and initialize with \a init.
334  */
335  void reshape(difference_type const & shape, const_reference init)
336  {
337  BaseType::reshape(shape, init);
338  }
339 
340  /** Create a matrix view that represents the row vector of row \a d.
341  */
343  {
344  return vigra::linalg::rowVector(*this, d);
345  }
346 
347  /** Create a matrix view that represents the column vector of column \a d.
348  */
350  {
351  return vigra::linalg::columnVector(*this, d);
352  }
353 
354  /** number of rows (height) of the matrix.
355  */
357  {
358  return this->m_shape[0];
359  }
360 
361  /** number of columns (width) of the matrix.
362  */
364  {
365  return this->m_shape[1];
366  }
367 
368  /** number of elements (width*height) of the matrix.
369  */
371  {
372  return rowCount()*columnCount();
373  }
374 
375  /** check whether the matrix is symmetric.
376  */
377  bool isSymmetric() const
378  {
379  return vigra::linalg::isSymmetric(*this);
380  }
381 
382 #ifdef DOXYGEN
383 // repeat the following functions for documentation. In real code, they are inherited.
384 
385  /** read/write access to matrix element <tt>(row, column)</tt>.
386  Note that the order of the argument is the opposite of the usual
387  VIGRA convention due to column-major matrix order.
388  */
390 
391  /** read access to matrix element <tt>(row, column)</tt>.
392  Note that the order of the argument is the opposite of the usual
393  VIGRA convention due to column-major matrix order.
394  */
396 
397  /** squared Frobenius norm. Sum of squares of the matrix elements.
398  */
399  typename NormTraits<Matrix>::SquaredNormType squaredNorm() const;
400 
401  /** Frobenius norm. Root of sum of squares of the matrix elements.
402  */
403  typename NormTraits<Matrix>::NormType norm() const;
404 
405  /** create a transposed view of this matrix.
406  No data are copied. If you want to transpose this matrix permanently,
407  you have to assign the transposed view:
408 
409  \code
410  a = a.transpose();
411  \endcode
412  */
414 #endif
415 
416  /** add \a other to this (sizes must match).
417  */
418  template <class U, class C>
420  {
421  BaseType::operator+=(other);
422  return *this;
423  }
424 
425  /** subtract \a other from this (sizes must match).
426  */
427  template <class U, class C>
429  {
430  BaseType::operator-=(other);
431  return *this;
432  }
433 
434  /** multiply \a other element-wise with this matrix (sizes must match).
435  */
436  template <class U, class C>
438  {
439  BaseType::operator*=(other);
440  return *this;
441  }
442 
443  /** divide this matrix element-wise by \a other (sizes must match).
444  */
445  template <class U, class C>
447  {
448  BaseType::operator/=(other);
449  return *this;
450  }
451 
452  /** add \a other to each element of this matrix
453  */
454  Matrix & operator+=(T other)
455  {
456  BaseType::operator+=(other);
457  return *this;
458  }
459 
460  /** subtraxt \a other from each element of this matrix
461  */
462  Matrix & operator-=(T other)
463  {
464  BaseType::operator-=(other);
465  return *this;
466  }
467 
468  /** scalar multiply this with \a other
469  */
470  Matrix & operator*=(T other)
471  {
472  BaseType::operator*=(other);
473  return *this;
474  }
475 
476  /** scalar devide this by \a other
477  */
478  Matrix & operator/=(T other)
479  {
480  BaseType::operator/=(other);
481  return *this;
482  }
483 };
484 
485 // TemporaryMatrix is provided as an optimization: Functions returning a matrix can
486 // use TemporaryMatrix to make explicit that it was allocated as a temporary data structure.
487 // Functions receiving a TemporaryMatrix can thus often avoid to allocate new temporary
488 // memory.
489 template <class T, class ALLOC = std::allocator<T> >
490 class TemporaryMatrix
491 : public Matrix<T, ALLOC>
492 {
493  typedef Matrix<T, ALLOC> BaseType;
494  public:
495  typedef Matrix<T, ALLOC> matrix_type;
496  typedef TemporaryMatrix<T, ALLOC> temp_type;
497  typedef MultiArrayView<2, T, UnstridedArrayTag> view_type;
498  typedef typename BaseType::value_type value_type;
499  typedef typename BaseType::pointer pointer;
500  typedef typename BaseType::const_pointer const_pointer;
501  typedef typename BaseType::reference reference;
502  typedef typename BaseType::const_reference const_reference;
503  typedef typename BaseType::difference_type difference_type;
504  typedef typename BaseType::difference_type_1 difference_type_1;
505  typedef ALLOC allocator_type;
506 
507  TemporaryMatrix(difference_type const & shape)
508  : BaseType(shape, ALLOC())
509  {}
510 
511  TemporaryMatrix(difference_type const & shape, const_reference init)
512  : BaseType(shape, init, ALLOC())
513  {}
514 
515  TemporaryMatrix(difference_type_1 rows, difference_type_1 columns)
516  : BaseType(rows, columns, ALLOC())
517  {}
518 
519  TemporaryMatrix(difference_type_1 rows, difference_type_1 columns, const_reference init)
520  : BaseType(rows, columns, init, ALLOC())
521  {}
522 
523  template<class U, class C>
524  TemporaryMatrix(const MultiArrayView<2, U, C> &rhs)
525  : BaseType(rhs)
526  {}
527 
528  TemporaryMatrix(const TemporaryMatrix &rhs)
529  : BaseType()
530  {
531  this->swap(const_cast<TemporaryMatrix &>(rhs));
532  }
533 
534  template <class U>
535  TemporaryMatrix & init(const U & init)
536  {
537  BaseType::init(init);
538  return *this;
539  }
540 
541  template <class U, class C>
542  TemporaryMatrix & operator+=(MultiArrayView<2, U, C> const & other)
543  {
544  BaseType::operator+=(other);
545  return *this;
546  }
547 
548  template <class U, class C>
549  TemporaryMatrix & operator-=(MultiArrayView<2, U, C> const & other)
550  {
551  BaseType::operator-=(other);
552  return *this;
553  }
554 
555  template <class U, class C>
556  TemporaryMatrix & operator*=(MultiArrayView<2, U, C> const & other)
557  {
558  BaseType::operator*=(other);
559  return *this;
560  }
561 
562  template <class U, class C>
563  TemporaryMatrix & operator/=(MultiArrayView<2, U, C> const & other)
564  {
565  BaseType::operator/=(other);
566  return *this;
567  }
568 
569  TemporaryMatrix & operator+=(T other)
570  {
571  BaseType::operator+=(other);
572  return *this;
573  }
574 
575  TemporaryMatrix & operator-=(T other)
576  {
577  BaseType::operator-=(other);
578  return *this;
579  }
580 
581  TemporaryMatrix & operator*=(T other)
582  {
583  BaseType::operator*=(other);
584  return *this;
585  }
586 
587  TemporaryMatrix & operator/=(T other)
588  {
589  BaseType::operator/=(other);
590  return *this;
591  }
592  private:
593 
594  TemporaryMatrix &operator=(const TemporaryMatrix &rhs); // intentionally not implemented
595 };
596 
597 /** \defgroup LinearAlgebraFunctions Matrix Functions
598 
599  \brief Basic matrix algebra, element-wise mathematical functions, row and columns statistics, data normalization etc.
600 
601  \ingroup LinearAlgebraModule
602  */
603 //@{
604 
605  /** Number of rows of a matrix represented as a <tt>MultiArrayView<2, ...></tt>
606 
607  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
608  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
609  Namespaces: vigra and vigra::linalg
610  */
611 template <class T, class C>
612 inline MultiArrayIndex
614 {
615  return x.shape(0);
616 }
617 
618  /** Number of columns of a matrix represented as a <tt>MultiArrayView<2, ...></tt>
619 
620  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
621  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
622  Namespaces: vigra and vigra::linalg
623  */
624 template <class T, class C>
625 inline MultiArrayIndex
627 {
628  return x.shape(1);
629 }
630 
631  /** Create a row vector view for row \a d of the matrix \a m
632 
633  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
634  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
635  Namespaces: vigra and vigra::linalg
636  */
637 template <class T, class C>
640 {
641  typedef typename MultiArrayView <2, T, C>::difference_type Shape;
642  return m.subarray(Shape(d, 0), Shape(d+1, columnCount(m)));
643 }
644 
645 
646  /** Create a row vector view of the matrix \a m starting at element \a first and ranging
647  to column \a end (non-inclusive).
648 
649  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
650  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
651  Namespaces: vigra and vigra::linalg
652  */
653 template <class T, class C>
656 {
657  typedef typename MultiArrayView <2, T, C>::difference_type Shape;
658  return m.subarray(first, Shape(first[0]+1, end));
659 }
660 
661  /** Create a column vector view for column \a d of the matrix \a m
662 
663  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
664  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
665  Namespaces: vigra and vigra::linalg
666  */
667 template <class T, class C>
670 {
671  typedef typename MultiArrayView <2, T, C>::difference_type Shape;
672  return m.subarray(Shape(0, d), Shape(rowCount(m), d+1));
673 }
674 
675  /** Create a column vector view of the matrix \a m starting at element \a first and
676  ranging to row \a end (non-inclusive).
677 
678  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
679  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
680  Namespaces: vigra and vigra::linalg
681  **/
682 template <class T, class C>
685 {
686  typedef typename MultiArrayView <2, T, C>::difference_type Shape;
687  return m.subarray(first, Shape(end, first[1]+1));
688 }
689 
690  /** Create a sub vector view of the vector \a m starting at element \a first and
691  ranging to row \a end (non-inclusive).
692 
693  Note: This function may only be called when either <tt>rowCount(m) == 1</tt> or
694  <tt>columnCount(m) == 1</tt>, i.e. when \a m really represents a vector.
695  Otherwise, a PreconditionViolation exception is raised.
696 
697  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
698  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
699  Namespaces: vigra and vigra::linalg
700  **/
701 template <class T, class C>
703 subVector(MultiArrayView<2, T, C> const & m, int first, int end)
704 {
705  typedef typename MultiArrayView <2, T, C>::difference_type Shape;
706  if(columnCount(m) == 1)
707  return m.subarray(Shape(first, 0), Shape(end, 1));
708  vigra_precondition(rowCount(m) == 1,
709  "linalg::subVector(): Input must be a vector (1xN or Nx1).");
710  return m.subarray(Shape(0, first), Shape(1, end));
711 }
712 
713  /** Check whether matrix \a m is symmetric.
714 
715  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
716  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
717  Namespaces: vigra and vigra::linalg
718  */
719 template <class T, class C>
720 bool
722 {
723  const MultiArrayIndex size = rowCount(m);
724  if(size != columnCount(m))
725  return false;
726 
727  for(MultiArrayIndex i = 0; i < size; ++i)
728  for(MultiArrayIndex j = i+1; j < size; ++j)
729  if(m(j, i) != m(i, j))
730  return false;
731  return true;
732 }
733 
734 
735  /** Compute the trace of a square matrix.
736 
737  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
738  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
739  Namespaces: vigra and vigra::linalg
740  */
741 template <class T, class C>
742 typename NumericTraits<T>::Promote
744 {
745  typedef typename NumericTraits<T>::Promote SumType;
746 
747  const MultiArrayIndex size = rowCount(m);
748  vigra_precondition(size == columnCount(m), "linalg::trace(): Matrix must be square.");
749 
750  SumType sum = NumericTraits<SumType>::zero();
751  for(MultiArrayIndex i = 0; i < size; ++i)
752  sum += m(i, i);
753  return sum;
754 }
755 
756 #ifdef DOXYGEN // documentation only -- function is already defined in vigra/multi_array.hxx
757 
758  /** calculate the squared Frobenius norm of a matrix.
759  Equal to the sum of squares of the matrix elements.
760 
761  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>>
762  Namespace: vigra
763  */
764 template <class T, class ALLOC>
765 typename Matrix<T, ALLLOC>::SquaredNormType
766 squaredNorm(const Matrix<T, ALLLOC> &a);
767 
768  /** calculate the Frobenius norm of a matrix.
769  Equal to the root of the sum of squares of the matrix elements.
770 
771  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>>
772  Namespace: vigra
773  */
774 template <class T, class ALLOC>
775 typename Matrix<T, ALLLOC>::NormType
776 norm(const Matrix<T, ALLLOC> &a);
777 
778 #endif // DOXYGEN
779 
780  /** initialize the given square matrix as an identity matrix.
781 
782  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
783  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
784  Namespaces: vigra and vigra::linalg
785  */
786 template <class T, class C>
788 {
789  const MultiArrayIndex rows = rowCount(r);
790  vigra_precondition(rows == columnCount(r),
791  "identityMatrix(): Matrix must be square.");
792  for(MultiArrayIndex i = 0; i < rows; ++i) {
793  for(MultiArrayIndex j = 0; j < rows; ++j)
794  r(j, i) = NumericTraits<T>::zero();
795  r(i, i) = NumericTraits<T>::one();
796  }
797 }
798 
799  /** create n identity matrix of the given size.
800  Usage:
801 
802  \code
803  vigra::Matrix<double> m = vigra::identityMatrix<double>(size);
804  \endcode
805 
806  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
807  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
808  Namespaces: vigra and vigra::linalg
809  */
810 template <class T>
811 TemporaryMatrix<T> identityMatrix(MultiArrayIndex size)
812 {
813  TemporaryMatrix<T> ret(size, size, NumericTraits<T>::zero());
814  for(MultiArrayIndex i = 0; i < size; ++i)
815  ret(i, i) = NumericTraits<T>::one();
816  return ret;
817 }
818 
819 template <class T, class C1, class C2>
820 void diagonalMatrixImpl(MultiArrayView<1, T, C1> const & v, MultiArrayView<2, T, C2> &r)
821 {
822  const MultiArrayIndex size = v.elementCount();
823  vigra_precondition(rowCount(r) == size && columnCount(r) == size,
824  "diagonalMatrix(): result must be a square matrix.");
825  for(MultiArrayIndex i = 0; i < size; ++i)
826  r(i, i) = v(i);
827 }
828 
829  /** make a diagonal matrix from a vector.
830  The vector is given as matrix \a v, which must either have a single
831  row or column. The result is witten into the square matrix \a r.
832 
833  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
834  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
835  Namespaces: vigra and vigra::linalg
836  */
837 template <class T, class C1, class C2>
839 {
840  vigra_precondition(rowCount(v) == 1 || columnCount(v) == 1,
841  "diagonalMatrix(): input must be a vector.");
842  r.init(NumericTraits<T>::zero());
843  if(rowCount(v) == 1)
844  diagonalMatrixImpl(v.bindInner(0), r);
845  else
846  diagonalMatrixImpl(v.bindOuter(0), r);
847 }
848 
849  /** create a diagonal matrix from a vector.
850  The vector is given as matrix \a v, which must either have a single
851  row or column. The result is returned as a temporary matrix.
852  Usage:
853 
854  \code
855  vigra::Matrix<double> v(1, len);
856  v = ...;
857 
858  vigra::Matrix<double> m = diagonalMatrix(v);
859  \endcode
860 
861  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
862  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
863  Namespaces: vigra and vigra::linalg
864  */
865 template <class T, class C>
866 TemporaryMatrix<T> diagonalMatrix(MultiArrayView<2, T, C> const & v)
867 {
868  vigra_precondition(rowCount(v) == 1 || columnCount(v) == 1,
869  "diagonalMatrix(): input must be a vector.");
870  MultiArrayIndex size = v.elementCount();
871  TemporaryMatrix<T> ret(size, size, NumericTraits<T>::zero());
872  if(rowCount(v) == 1)
873  diagonalMatrixImpl(v.bindInner(0), ret);
874  else
875  diagonalMatrixImpl(v.bindOuter(0), ret);
876  return ret;
877 }
878 
879  /** transpose matrix \a v.
880  The result is written into \a r which must have the correct (i.e.
881  transposed) shape.
882 
883  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
884  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
885  Namespaces: vigra and vigra::linalg
886  */
887 template <class T, class C1, class C2>
889 {
890  const MultiArrayIndex rows = rowCount(r);
891  const MultiArrayIndex cols = columnCount(r);
892  vigra_precondition(rows == columnCount(v) && cols == rowCount(v),
893  "transpose(): arrays must have transposed shapes.");
894  for(MultiArrayIndex i = 0; i < cols; ++i)
895  for(MultiArrayIndex j = 0; j < rows; ++j)
896  r(j, i) = v(i, j);
897 }
898 
899  /** create the transpose of matrix \a v.
900  This does not copy any data, but only creates a transposed view
901  to the original matrix. A copy is only made when the transposed view
902  is assigned to another matrix.
903  Usage:
904 
905  \code
906  vigra::Matrix<double> v(rows, cols);
907  v = ...;
908 
909  vigra::Matrix<double> m = transpose(v);
910  \endcode
911 
912  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
913  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
914  Namespaces: vigra and vigra::linalg
915  */
916 template <class T, class C>
919 {
920  return v.transpose();
921 }
922 
923  /** Create new matrix by concatenating two matrices \a a and \a b vertically, i.e. on top of each other.
924  The two matrices must have the same number of columns.
925  The result is returned as a temporary matrix.
926 
927  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
928  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
929  Namespace: vigra::linalg
930  */
931 template <class T, class C1, class C2>
932 inline TemporaryMatrix<T>
934 {
935  typedef typename TemporaryMatrix<T>::difference_type Shape;
936 
938  vigra_precondition(n == columnCount(b),
939  "joinVertically(): shape mismatch.");
940 
941  MultiArrayIndex ma = rowCount(a);
942  MultiArrayIndex mb = rowCount(b);
943  TemporaryMatrix<T> t(ma + mb, n, T());
944  t.subarray(Shape(0,0), Shape(ma, n)) = a;
945  t.subarray(Shape(ma,0), Shape(ma+mb, n)) = b;
946  return t;
947 }
948 
949  /** Create new matrix by concatenating two matrices \a a and \a b horizontally, i.e. side by side.
950  The two matrices must have the same number of rows.
951  The result is returned as a temporary matrix.
952 
953  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
954  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
955  Namespace: vigra::linalg
956  */
957 template <class T, class C1, class C2>
958 inline TemporaryMatrix<T>
960 {
961  typedef typename TemporaryMatrix<T>::difference_type Shape;
962 
963  MultiArrayIndex m = rowCount(a);
964  vigra_precondition(m == rowCount(b),
965  "joinHorizontally(): shape mismatch.");
966 
969  TemporaryMatrix<T> t(m, na + nb, T());
970  t.subarray(Shape(0,0), Shape(m, na)) = a;
971  t.subarray(Shape(0, na), Shape(m, na + nb)) = b;
972  return t;
973 }
974 
975  /** Initialize a matrix with repeated copies of a given matrix.
976 
977  Matrix \a r will consist of \a verticalCount downward repetitions of \a v,
978  and \a horizontalCount side-by-side repetitions. When \a v has size <tt>m</tt> by <tt>n</tt>,
979  \a r must have size <tt>(m*verticalCount)</tt> by <tt>(n*horizontalCount)</tt>.
980 
981  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
982  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
983  Namespace: vigra::linalg
984  */
985 template <class T, class C1, class C2>
987  unsigned int verticalCount, unsigned int horizontalCount)
988 {
989  typedef typename Matrix<T>::difference_type Shape;
990 
991  MultiArrayIndex m = rowCount(v), n = columnCount(v);
992  vigra_precondition(m*verticalCount == rowCount(r) && n*horizontalCount == columnCount(r),
993  "repeatMatrix(): Shape mismatch.");
994 
995  for(MultiArrayIndex l=0; l<(MultiArrayIndex)horizontalCount; ++l)
996  {
997  for(MultiArrayIndex k=0; k<(MultiArrayIndex)verticalCount; ++k)
998  {
999  r.subarray(Shape(k*m, l*n), Shape((k+1)*m, (l+1)*n)) = v;
1000  }
1001  }
1002 }
1003 
1004  /** Create a new matrix by repeating a given matrix.
1005 
1006  The resulting matrix \a r will consist of \a verticalCount downward repetitions of \a v,
1007  and \a horizontalCount side-by-side repetitions, i.e. it will be of size
1008  <tt>(m*verticalCount)</tt> by <tt>(n*horizontalCount)</tt> when \a v has size <tt>m</tt> by <tt>n</tt>.
1009  The result is returned as a temporary matrix.
1010 
1011  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1012  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1013  Namespace: vigra::linalg
1014  */
1015 template <class T, class C>
1016 TemporaryMatrix<T>
1017 repeatMatrix(MultiArrayView<2, T, C> const & v, unsigned int verticalCount, unsigned int horizontalCount)
1018 {
1019  MultiArrayIndex m = rowCount(v), n = columnCount(v);
1020  TemporaryMatrix<T> ret(verticalCount*m, horizontalCount*n);
1021  repeatMatrix(v, ret, verticalCount, horizontalCount);
1022  return ret;
1023 }
1024 
1025  /** add matrices \a a and \a b.
1026  The result is written into \a r. All three matrices must have the same shape.
1027 
1028  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1029  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1030  Namespace: vigra::linalg
1031  */
1032 template <class T, class C1, class C2, class C3>
1035 {
1036  const MultiArrayIndex rrows = rowCount(r);
1037  const MultiArrayIndex rcols = columnCount(r);
1038  vigra_precondition(rrows == rowCount(a) && rcols == columnCount(a) &&
1039  rrows == rowCount(b) && rcols == columnCount(b),
1040  "add(): Matrix shapes must agree.");
1041 
1042  for(MultiArrayIndex i = 0; i < rcols; ++i) {
1043  for(MultiArrayIndex j = 0; j < rrows; ++j) {
1044  r(j, i) = a(j, i) + b(j, i);
1045  }
1046  }
1047 }
1048 
1049  /** add matrices \a a and \a b.
1050  The two matrices must have the same shape.
1051  The result is returned as a temporary matrix.
1052 
1053  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1054  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1055  Namespace: vigra::linalg
1056  */
1057 template <class T, class C1, class C2>
1058 inline TemporaryMatrix<T>
1060 {
1061  return TemporaryMatrix<T>(a) += b;
1062 }
1063 
1064 template <class T, class C>
1065 inline TemporaryMatrix<T>
1066 operator+(const TemporaryMatrix<T> &a, const MultiArrayView<2, T, C> &b)
1067 {
1068  return const_cast<TemporaryMatrix<T> &>(a) += b;
1069 }
1070 
1071 template <class T, class C>
1072 inline TemporaryMatrix<T>
1073 operator+(const MultiArrayView<2, T, C> &a, const TemporaryMatrix<T> &b)
1074 {
1075  return const_cast<TemporaryMatrix<T> &>(b) += a;
1076 }
1077 
1078 template <class T>
1079 inline TemporaryMatrix<T>
1080 operator+(const TemporaryMatrix<T> &a, const TemporaryMatrix<T> &b)
1081 {
1082  return const_cast<TemporaryMatrix<T> &>(a) += b;
1083 }
1084 
1085  /** add scalar \a b to every element of the matrix \a a.
1086  The result is returned as a temporary matrix.
1087 
1088  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1089  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1090  Namespace: vigra::linalg
1091  */
1092 template <class T, class C>
1093 inline TemporaryMatrix<T>
1095 {
1096  return TemporaryMatrix<T>(a) += b;
1097 }
1098 
1099 template <class T>
1100 inline TemporaryMatrix<T>
1101 operator+(const TemporaryMatrix<T> &a, T b)
1102 {
1103  return const_cast<TemporaryMatrix<T> &>(a) += b;
1104 }
1105 
1106  /** add scalar \a a to every element of the matrix \a b.
1107  The result is returned as a temporary matrix.
1108 
1109  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1110  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1111  Namespace: vigra::linalg
1112  */
1113 template <class T, class C>
1114 inline TemporaryMatrix<T>
1116 {
1117  return TemporaryMatrix<T>(b) += a;
1118 }
1119 
1120 template <class T>
1121 inline TemporaryMatrix<T>
1122 operator+(T a, const TemporaryMatrix<T> &b)
1123 {
1124  return const_cast<TemporaryMatrix<T> &>(b) += a;
1125 }
1126 
1127  /** subtract matrix \a b from \a a.
1128  The result is written into \a r. All three matrices must have the same shape.
1129 
1130  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1131  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1132  Namespace: vigra::linalg
1133  */
1134 template <class T, class C1, class C2, class C3>
1137 {
1138  const MultiArrayIndex rrows = rowCount(r);
1139  const MultiArrayIndex rcols = columnCount(r);
1140  vigra_precondition(rrows == rowCount(a) && rcols == columnCount(a) &&
1141  rrows == rowCount(b) && rcols == columnCount(b),
1142  "subtract(): Matrix shapes must agree.");
1143 
1144  for(MultiArrayIndex i = 0; i < rcols; ++i) {
1145  for(MultiArrayIndex j = 0; j < rrows; ++j) {
1146  r(j, i) = a(j, i) - b(j, i);
1147  }
1148  }
1149 }
1150 
1151  /** subtract matrix \a b from \a a.
1152  The two matrices must have the same shape.
1153  The result is returned as a temporary matrix.
1154 
1155  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1156  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1157  Namespace: vigra::linalg
1158  */
1159 template <class T, class C1, class C2>
1160 inline TemporaryMatrix<T>
1162 {
1163  return TemporaryMatrix<T>(a) -= b;
1164 }
1165 
1166 template <class T, class C>
1167 inline TemporaryMatrix<T>
1168 operator-(const TemporaryMatrix<T> &a, const MultiArrayView<2, T, C> &b)
1169 {
1170  return const_cast<TemporaryMatrix<T> &>(a) -= b;
1171 }
1172 
1173 template <class T, class C>
1174 TemporaryMatrix<T>
1175 operator-(const MultiArrayView<2, T, C> &a, const TemporaryMatrix<T> &b)
1176 {
1177  const MultiArrayIndex rows = rowCount(a);
1178  const MultiArrayIndex cols = columnCount(a);
1179  vigra_precondition(rows == b.rowCount() && cols == b.columnCount(),
1180  "Matrix::operator-(): Shape mismatch.");
1181 
1182  for(MultiArrayIndex i = 0; i < cols; ++i)
1183  for(MultiArrayIndex j = 0; j < rows; ++j)
1184  const_cast<TemporaryMatrix<T> &>(b)(j, i) = a(j, i) - b(j, i);
1185  return b;
1186 }
1187 
1188 template <class T>
1189 inline TemporaryMatrix<T>
1190 operator-(const TemporaryMatrix<T> &a, const TemporaryMatrix<T> &b)
1191 {
1192  return const_cast<TemporaryMatrix<T> &>(a) -= b;
1193 }
1194 
1195  /** negate matrix \a a.
1196  The result is returned as a temporary matrix.
1197 
1198  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1199  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1200  Namespace: vigra::linalg
1201  */
1202 template <class T, class C>
1203 inline TemporaryMatrix<T>
1205 {
1206  return TemporaryMatrix<T>(a) *= -NumericTraits<T>::one();
1207 }
1208 
1209 template <class T>
1210 inline TemporaryMatrix<T>
1211 operator-(const TemporaryMatrix<T> &a)
1212 {
1213  return const_cast<TemporaryMatrix<T> &>(a) *= -NumericTraits<T>::one();
1214 }
1215 
1216  /** subtract scalar \a b from every element of the matrix \a a.
1217  The result is returned as a temporary matrix.
1218 
1219  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1220  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1221  Namespace: vigra::linalg
1222  */
1223 template <class T, class C>
1224 inline TemporaryMatrix<T>
1226 {
1227  return TemporaryMatrix<T>(a) -= b;
1228 }
1229 
1230 template <class T>
1231 inline TemporaryMatrix<T>
1232 operator-(const TemporaryMatrix<T> &a, T b)
1233 {
1234  return const_cast<TemporaryMatrix<T> &>(a) -= b;
1235 }
1236 
1237  /** subtract every element of the matrix \a b from scalar \a a.
1238  The result is returned as a temporary matrix.
1239 
1240  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1241  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1242  Namespace: vigra::linalg
1243  */
1244 template <class T, class C>
1245 inline TemporaryMatrix<T>
1247 {
1248  return TemporaryMatrix<T>(b.shape(), a) -= b;
1249 }
1250 
1251  /** calculate the inner product of two matrices representing vectors.
1252  Typically, matrix \a x has a single row, and matrix \a y has
1253  a single column, and the other dimensions match. In addition, this
1254  function handles the cases when either or both of the two inputs are
1255  transposed (e.g. it can compute the dot product of two column vectors).
1256  A <tt>PreconditionViolation</tt> exception is thrown when
1257  the shape conditions are violated.
1258 
1259  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1260  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1261  Namespaces: vigra and vigra::linalg
1262  */
1263 template <class T, class C1, class C2>
1264 typename NormTraits<T>::SquaredNormType
1266 {
1267  typename NormTraits<T>::SquaredNormType ret =
1268  NumericTraits<typename NormTraits<T>::SquaredNormType>::zero();
1269  if(y.shape(1) == 1)
1270  {
1271  std::ptrdiff_t size = y.shape(0);
1272  if(x.shape(0) == 1 && x.shape(1) == size) // proper scalar product
1273  for(std::ptrdiff_t i = 0; i < size; ++i)
1274  ret += x(0, i) * y(i, 0);
1275  else if(x.shape(1) == 1u && x.shape(0) == size) // two column vectors
1276  for(std::ptrdiff_t i = 0; i < size; ++i)
1277  ret += x(i, 0) * y(i, 0);
1278  else
1279  vigra_precondition(false, "dot(): wrong matrix shapes.");
1280  }
1281  else if(y.shape(0) == 1)
1282  {
1283  std::ptrdiff_t size = y.shape(1);
1284  if(x.shape(0) == 1u && x.shape(1) == size) // two row vectors
1285  for(std::ptrdiff_t i = 0; i < size; ++i)
1286  ret += x(0, i) * y(0, i);
1287  else if(x.shape(1) == 1u && x.shape(0) == size) // column dot row
1288  for(std::ptrdiff_t i = 0; i < size; ++i)
1289  ret += x(i, 0) * y(0, i);
1290  else
1291  vigra_precondition(false, "dot(): wrong matrix shapes.");
1292  }
1293  else
1294  vigra_precondition(false, "dot(): wrong matrix shapes.");
1295  return ret;
1296 }
1297 
1298  /** calculate the inner product of two vectors. The vector
1299  lengths must match.
1300 
1301  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1302  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1303  Namespaces: vigra and vigra::linalg
1304  */
1305 template <class T, class C1, class C2>
1306 typename NormTraits<T>::SquaredNormType
1308 {
1309  const MultiArrayIndex n = x.elementCount();
1310  vigra_precondition(n == y.elementCount(),
1311  "dot(): shape mismatch.");
1312  typename NormTraits<T>::SquaredNormType ret =
1313  NumericTraits<typename NormTraits<T>::SquaredNormType>::zero();
1314  for(MultiArrayIndex i = 0; i < n; ++i)
1315  ret += x(i) * y(i);
1316  return ret;
1317 }
1318 
1319  /** calculate the cross product of two vectors of length 3.
1320  The result is written into \a r.
1321 
1322  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1323  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1324  Namespaces: vigra and vigra::linalg
1325  */
1326 template <class T, class C1, class C2, class C3>
1329 {
1330  vigra_precondition(3 == x.elementCount() && 3 == y.elementCount() && 3 == r.elementCount(),
1331  "cross(): vectors must have length 3.");
1332  r(0) = x(1)*y(2) - x(2)*y(1);
1333  r(1) = x(2)*y(0) - x(0)*y(2);
1334  r(2) = x(0)*y(1) - x(1)*y(0);
1335 }
1336 
1337  /** calculate the cross product of two matrices representing vectors.
1338  That is, \a x, \a y, and \a r must have a single column of length 3. The result
1339  is written into \a r.
1340 
1341  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1342  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1343  Namespaces: vigra and vigra::linalg
1344  */
1345 template <class T, class C1, class C2, class C3>
1348 {
1349  vigra_precondition(3 == rowCount(x) && 3 == rowCount(y) && 3 == rowCount(r) ,
1350  "cross(): vectors must have length 3.");
1351  r(0,0) = x(1,0)*y(2,0) - x(2,0)*y(1,0);
1352  r(1,0) = x(2,0)*y(0,0) - x(0,0)*y(2,0);
1353  r(2,0) = x(0,0)*y(1,0) - x(1,0)*y(0,0);
1354 }
1355 
1356  /** calculate the cross product of two matrices representing vectors.
1357  That is, \a x, and \a y must have a single column of length 3. The result
1358  is returned as a temporary matrix.
1359 
1360  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1361  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1362  Namespaces: vigra and vigra::linalg
1363  */
1364 template <class T, class C1, class C2>
1365 TemporaryMatrix<T>
1367 {
1368  TemporaryMatrix<T> ret(3, 1);
1369  cross(x, y, ret);
1370  return ret;
1371 }
1372  /** calculate the outer product of two matrices representing vectors.
1373  That is, matrix \a x must have a single column, and matrix \a y must
1374  have a single row, and the other dimensions must match. The result
1375  is written into \a r.
1376 
1377  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1378  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1379  Namespaces: vigra and vigra::linalg
1380  */
1381 template <class T, class C1, class C2, class C3>
1384 {
1385  const MultiArrayIndex rows = rowCount(r);
1386  const MultiArrayIndex cols = columnCount(r);
1387  vigra_precondition(rows == rowCount(x) && cols == columnCount(y) &&
1388  1 == columnCount(x) && 1 == rowCount(y),
1389  "outer(): shape mismatch.");
1390  for(MultiArrayIndex i = 0; i < cols; ++i)
1391  for(MultiArrayIndex j = 0; j < rows; ++j)
1392  r(j, i) = x(j, 0) * y(0, i);
1393 }
1394 
1395  /** calculate the outer product of two matrices representing vectors.
1396  That is, matrix \a x must have a single column, and matrix \a y must
1397  have a single row, and the other dimensions must match. The result
1398  is returned as a temporary matrix.
1399 
1400  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1401  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1402  Namespaces: vigra and vigra::linalg
1403  */
1404 template <class T, class C1, class C2>
1405 TemporaryMatrix<T>
1407 {
1408  const MultiArrayIndex rows = rowCount(x);
1409  const MultiArrayIndex cols = columnCount(y);
1410  vigra_precondition(1 == columnCount(x) && 1 == rowCount(y),
1411  "outer(): shape mismatch.");
1412  TemporaryMatrix<T> ret(rows, cols);
1413  outer(x, y, ret);
1414  return ret;
1415 }
1416 
1417  /** calculate the outer product of a matrix (representing a vector) with itself.
1418  The result is returned as a temporary matrix.
1419 
1420  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1421  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1422  Namespaces: vigra and vigra::linalg
1423  */
1424 template <class T, class C>
1425 TemporaryMatrix<T>
1427 {
1428  const MultiArrayIndex rows = rowCount(x);
1429  const MultiArrayIndex cols = columnCount(x);
1430  vigra_precondition(rows == 1 || cols == 1,
1431  "outer(): matrix does not represent a vector.");
1432  const MultiArrayIndex size = std::max(rows, cols);
1433  TemporaryMatrix<T> ret(size, size);
1434 
1435  if(rows == 1)
1436  {
1437  for(MultiArrayIndex i = 0; i < size; ++i)
1438  for(MultiArrayIndex j = 0; j < size; ++j)
1439  ret(j, i) = x(0, j) * x(0, i);
1440  }
1441  else
1442  {
1443  for(MultiArrayIndex i = 0; i < size; ++i)
1444  for(MultiArrayIndex j = 0; j < size; ++j)
1445  ret(j, i) = x(j, 0) * x(i, 0);
1446  }
1447  return ret;
1448 }
1449 
1450 template <class T>
1451 class PointWise
1452 {
1453  public:
1454  T const & t;
1455 
1456  PointWise(T const & it)
1457  : t(it)
1458  {}
1459 };
1460 
1461 template <class T>
1462 PointWise<T> pointWise(T const & t)
1463 {
1464  return PointWise<T>(t);
1465 }
1466 
1467 
1468  /** multiply matrix \a a with scalar \a b.
1469  The result is written into \a r. \a a and \a r must have the same shape.
1470 
1471  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1472  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1473  Namespace: vigra::linalg
1474  */
1475 template <class T, class C1, class C2>
1477 {
1478  const MultiArrayIndex rows = rowCount(a);
1479  const MultiArrayIndex cols = columnCount(a);
1480  vigra_precondition(rows == rowCount(r) && cols == columnCount(r),
1481  "smul(): Matrix sizes must agree.");
1482 
1483  for(MultiArrayIndex i = 0; i < cols; ++i)
1484  for(MultiArrayIndex j = 0; j < rows; ++j)
1485  r(j, i) = a(j, i) * b;
1486 }
1487 
1488  /** multiply scalar \a a with matrix \a b.
1489  The result is written into \a r. \a b and \a r must have the same shape.
1490 
1491  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1492  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1493  Namespace: vigra::linalg
1494  */
1495 template <class T, class C2, class C3>
1497 {
1498  smul(b, a, r);
1499 }
1500 
1501  /** perform matrix multiplication of matrices \a a and \a b.
1502  The result is written into \a r. The three matrices must have matching shapes.
1503 
1504  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1505  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1506  Namespace: vigra::linalg
1507  */
1508 template <class T, class C1, class C2, class C3>
1511 {
1512  const MultiArrayIndex rrows = rowCount(r);
1513  const MultiArrayIndex rcols = columnCount(r);
1514  const MultiArrayIndex acols = columnCount(a);
1515  vigra_precondition(rrows == rowCount(a) && rcols == columnCount(b) && acols == rowCount(b),
1516  "mmul(): Matrix shapes must agree.");
1517 
1518  // order of loops ensures that inner loop goes down columns
1519  for(MultiArrayIndex i = 0; i < rcols; ++i)
1520  {
1521  for(MultiArrayIndex j = 0; j < rrows; ++j)
1522  r(j, i) = a(j, 0) * b(0, i);
1523  for(MultiArrayIndex k = 1; k < acols; ++k)
1524  for(MultiArrayIndex j = 0; j < rrows; ++j)
1525  r(j, i) += a(j, k) * b(k, i);
1526  }
1527 }
1528 
1529  /** perform matrix multiplication of matrices \a a and \a b.
1530  \a a and \a b must have matching shapes.
1531  The result is returned as a temporary matrix.
1532 
1533  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1534  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1535  Namespace: vigra::linalg
1536  */
1537 template <class T, class C1, class C2>
1538 inline TemporaryMatrix<T>
1540 {
1541  TemporaryMatrix<T> ret(rowCount(a), columnCount(b));
1542  mmul(a, b, ret);
1543  return ret;
1544 }
1545 
1546  /** multiply two matrices \a a and \a b pointwise.
1547  The result is written into \a r. All three matrices must have the same shape.
1548 
1549  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1550  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1551  Namespace: vigra::linalg
1552  */
1553 template <class T, class C1, class C2, class C3>
1556 {
1557  const MultiArrayIndex rrows = rowCount(r);
1558  const MultiArrayIndex rcols = columnCount(r);
1559  vigra_precondition(rrows == rowCount(a) && rcols == columnCount(a) &&
1560  rrows == rowCount(b) && rcols == columnCount(b),
1561  "pmul(): Matrix shapes must agree.");
1562 
1563  for(MultiArrayIndex i = 0; i < rcols; ++i) {
1564  for(MultiArrayIndex j = 0; j < rrows; ++j) {
1565  r(j, i) = a(j, i) * b(j, i);
1566  }
1567  }
1568 }
1569 
1570  /** multiply matrices \a a and \a b pointwise.
1571  \a a and \a b must have matching shapes.
1572  The result is returned as a temporary matrix.
1573 
1574  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1575  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1576  Namespace: vigra::linalg
1577  */
1578 template <class T, class C1, class C2>
1579 inline TemporaryMatrix<T>
1581 {
1582  TemporaryMatrix<T> ret(a.shape());
1583  pmul(a, b, ret);
1584  return ret;
1585 }
1586 
1587  /** multiply matrices \a a and \a b pointwise.
1588  \a a and \a b must have matching shapes.
1589  The result is returned as a temporary matrix.
1590 
1591  Usage:
1592 
1593  \code
1594  Matrix<double> a(m,n), b(m,n);
1595 
1596  Matrix<double> c = a * pointWise(b);
1597  // is equivalent to
1598  // Matrix<double> c = pmul(a, b);
1599  \endcode
1600 
1601  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1602  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1603  Namespace: vigra::linalg
1604  */
1605 template <class T, class C, class U>
1606 inline TemporaryMatrix<T>
1607 operator*(const MultiArrayView<2, T, C> &a, PointWise<U> b)
1608 {
1609  return pmul(a, b.t);
1610 }
1611 
1612  /** multiply matrix \a a with scalar \a b.
1613  The result is returned as a temporary matrix.
1614 
1615  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1616  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1617  Namespace: vigra::linalg
1618  */
1619 template <class T, class C>
1620 inline TemporaryMatrix<T>
1622 {
1623  return TemporaryMatrix<T>(a) *= b;
1624 }
1625 
1626 template <class T>
1627 inline TemporaryMatrix<T>
1628 operator*(const TemporaryMatrix<T> &a, T b)
1629 {
1630  return const_cast<TemporaryMatrix<T> &>(a) *= b;
1631 }
1632 
1633  /** multiply scalar \a a with matrix \a b.
1634  The result is returned as a temporary matrix.
1635 
1636  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1637  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1638  Namespace: vigra::linalg
1639  */
1640 template <class T, class C>
1641 inline TemporaryMatrix<T>
1643 {
1644  return TemporaryMatrix<T>(b) *= a;
1645 }
1646 
1647 template <class T>
1648 inline TemporaryMatrix<T>
1649 operator*(T a, const TemporaryMatrix<T> &b)
1650 {
1651  return const_cast<TemporaryMatrix<T> &>(b) *= a;
1652 }
1653 
1654  /** multiply matrix \a a with TinyVector \a b.
1655  \a a must be of size <tt>N x N</tt>. Vector \a b and the result
1656  vector are interpreted as column vectors.
1657 
1658  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1659  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1660  Namespace: vigra::linalg
1661  */
1662 template <class T, class A, int N, class DATA, class DERIVED>
1663 TinyVector<T, N>
1665 {
1666  vigra_precondition(N == rowCount(a) && N == columnCount(a),
1667  "operator*(Matrix, TinyVector): Shape mismatch.");
1668 
1669  TinyVector<T, N> res = TinyVectorView<T, N>(&a(0,0)) * b[0];
1670  for(MultiArrayIndex i = 1; i < N; ++i)
1671  res += TinyVectorView<T, N>(&a(0,i)) * b[i];
1672  return res;
1673 }
1674 
1675  /** multiply TinyVector \a a with matrix \a b.
1676  \a b must be of size <tt>N x N</tt>. Vector \a a and the result
1677  vector are interpreted as row vectors.
1678 
1679  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1680  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1681  Namespace: vigra::linalg
1682  */
1683 template <class T, int N, class DATA, class DERIVED, class A>
1686 {
1687  vigra_precondition(N == rowCount(b) && N == columnCount(b),
1688  "operator*(TinyVector, Matrix): Shape mismatch.");
1689 
1690  TinyVector<T, N> res;
1691  for(MultiArrayIndex i = 0; i < N; ++i)
1692  res[i] = dot(a, TinyVectorView<T, N>(&b(0,i)));
1693  return res;
1694 }
1695 
1696  /** perform matrix multiplication of matrices \a a and \a b.
1697  \a a and \a b must have matching shapes.
1698  The result is returned as a temporary matrix.
1699 
1700  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1701  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1702  Namespace: vigra::linalg
1703  */
1704 template <class T, class C1, class C2>
1705 inline TemporaryMatrix<T>
1707 {
1708  TemporaryMatrix<T> ret(rowCount(a), columnCount(b));
1709  mmul(a, b, ret);
1710  return ret;
1711 }
1712 
1713  /** divide matrix \a a by scalar \a b.
1714  The result is written into \a r. \a a and \a r must have the same shape.
1715 
1716  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1717  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1718  Namespace: vigra::linalg
1719  */
1720 template <class T, class C1, class C2>
1722 {
1723  const MultiArrayIndex rows = rowCount(a);
1724  const MultiArrayIndex cols = columnCount(a);
1725  vigra_precondition(rows == rowCount(r) && cols == columnCount(r),
1726  "sdiv(): Matrix sizes must agree.");
1727 
1728  for(MultiArrayIndex i = 0; i < cols; ++i)
1729  for(MultiArrayIndex j = 0; j < rows; ++j)
1730  r(j, i) = a(j, i) / b;
1731 }
1732 
1733  /** divide two matrices \a a and \a b pointwise.
1734  The result is written into \a r. All three matrices must have the same shape.
1735 
1736  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1737  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1738  Namespace: vigra::linalg
1739  */
1740 template <class T, class C1, class C2, class C3>
1743 {
1744  const MultiArrayIndex rrows = rowCount(r);
1745  const MultiArrayIndex rcols = columnCount(r);
1746  vigra_precondition(rrows == rowCount(a) && rcols == columnCount(a) &&
1747  rrows == rowCount(b) && rcols == columnCount(b),
1748  "pdiv(): Matrix shapes must agree.");
1749 
1750  for(MultiArrayIndex i = 0; i < rcols; ++i) {
1751  for(MultiArrayIndex j = 0; j < rrows; ++j) {
1752  r(j, i) = a(j, i) / b(j, i);
1753  }
1754  }
1755 }
1756 
1757  /** divide matrices \a a and \a b pointwise.
1758  \a a and \a b must have matching shapes.
1759  The result is returned as a temporary matrix.
1760 
1761  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1762  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1763  Namespace: vigra::linalg
1764  */
1765 template <class T, class C1, class C2>
1766 inline TemporaryMatrix<T>
1768 {
1769  TemporaryMatrix<T> ret(a.shape());
1770  pdiv(a, b, ret);
1771  return ret;
1772 }
1773 
1774  /** divide matrices \a a and \a b pointwise.
1775  \a a and \a b must have matching shapes.
1776  The result is returned as a temporary matrix.
1777 
1778  Usage:
1779 
1780  \code
1781  Matrix<double> a(m,n), b(m,n);
1782 
1783  Matrix<double> c = a / pointWise(b);
1784  // is equivalent to
1785  // Matrix<double> c = pdiv(a, b);
1786  \endcode
1787 
1788  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1789  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1790  Namespace: vigra::linalg
1791  */
1792 template <class T, class C, class U>
1793 inline TemporaryMatrix<T>
1794 operator/(const MultiArrayView<2, T, C> &a, PointWise<U> b)
1795 {
1796  return pdiv(a, b.t);
1797 }
1798 
1799  /** divide matrix \a a by scalar \a b.
1800  The result is returned as a temporary matrix.
1801 
1802  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1803  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1804  Namespace: vigra::linalg
1805  */
1806 template <class T, class C>
1807 inline TemporaryMatrix<T>
1809 {
1810  return TemporaryMatrix<T>(a) /= b;
1811 }
1812 
1813 template <class T>
1814 inline TemporaryMatrix<T>
1815 operator/(const TemporaryMatrix<T> &a, T b)
1816 {
1817  return const_cast<TemporaryMatrix<T> &>(a) /= b;
1818 }
1819 
1820  /** Create a matrix whose elements are the quotients between scalar \a a and
1821  matrix \a b. The result is returned as a temporary matrix.
1822 
1823  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
1824  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
1825  Namespace: vigra::linalg
1826  */
1827 template <class T, class C>
1828 inline TemporaryMatrix<T>
1830 {
1831  return TemporaryMatrix<T>(b.shape(), a) / pointWise(b);
1832 }
1833 
1834 using vigra::argMin;
1835 using vigra::argMinIf;
1836 using vigra::argMax;
1837 using vigra::argMaxIf;
1838 
1839  /*! Find the index of the minimum element in a matrix.
1840 
1841  The function returns the index in column-major scan-order sense,
1842  i.e. according to the order used by <tt>MultiArrayView::operator[]</tt>.
1843  If the matrix is actually a vector, this is just the row or columns index.
1844  In case of a truely 2-dimensional matrix, the index can be converted to an
1845  index pair by calling <tt>MultiArrayView::scanOrderIndexToCoordinate()</tt>
1846 
1847  <b>Required Interface:</b>
1848 
1849  \code
1850  bool f = a[0] < NumericTraits<T>::max();
1851  \endcode
1852 
1853  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>><br>
1854  Namespace: vigra
1855  */
1856 template <class T, class C>
1858 {
1859  T vopt = NumericTraits<T>::max();
1860  int best = -1;
1861  for(int k=0; k < a.size(); ++k)
1862  {
1863  if(a[k] < vopt)
1864  {
1865  vopt = a[k];
1866  best = k;
1867  }
1868  }
1869  return best;
1870 }
1871 
1872  /*! Find the index of the maximum element in a matrix.
1873 
1874  The function returns the index in column-major scan-order sense,
1875  i.e. according to the order used by <tt>MultiArrayView::operator[]</tt>.
1876  If the matrix is actually a vector, this is just the row or columns index.
1877  In case of a truely 2-dimensional matrix, the index can be converted to an
1878  index pair by calling <tt>MultiArrayView::scanOrderIndexToCoordinate()</tt>
1879 
1880  <b>Required Interface:</b>
1881 
1882  \code
1883  bool f = NumericTraits<T>::min() < a[0];
1884  \endcode
1885 
1886  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>><br>
1887  Namespace: vigra
1888  */
1889 template <class T, class C>
1891 {
1892  T vopt = NumericTraits<T>::min();
1893  int best = -1;
1894  for(int k=0; k < a.size(); ++k)
1895  {
1896  if(vopt < a[k])
1897  {
1898  vopt = a[k];
1899  best = k;
1900  }
1901  }
1902  return best;
1903 }
1904 
1905  /*! Find the index of the minimum element in a matrix subject to a condition.
1906 
1907  The function returns <tt>-1</tt> if no element conforms to \a condition.
1908  Otherwise, the index of the maximum element is returned in column-major scan-order sense,
1909  i.e. according to the order used by <tt>MultiArrayView::operator[]</tt>.
1910  If the matrix is actually a vector, this is just the row or columns index.
1911  In case of a truely 2-dimensional matrix, the index can be converted to an
1912  index pair by calling <tt>MultiArrayView::scanOrderIndexToCoordinate()</tt>
1913 
1914  <b>Required Interface:</b>
1915 
1916  \code
1917  bool c = condition(a[0]);
1918  bool f = a[0] < NumericTraits<T>::max();
1919  \endcode
1920 
1921  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>><br>
1922  Namespace: vigra
1923  */
1924 template <class T, class C, class UnaryFunctor>
1925 int argMinIf(MultiArrayView<2, T, C> const & a, UnaryFunctor condition)
1926 {
1927  T vopt = NumericTraits<T>::max();
1928  int best = -1;
1929  for(int k=0; k < a.size(); ++k)
1930  {
1931  if(condition(a[k]) && a[k] < vopt)
1932  {
1933  vopt = a[k];
1934  best = k;
1935  }
1936  }
1937  return best;
1938 }
1939 
1940  /*! Find the index of the maximum element in a matrix subject to a condition.
1941 
1942  The function returns <tt>-1</tt> if no element conforms to \a condition.
1943  Otherwise, the index of the maximum element is returned in column-major scan-order sense,
1944  i.e. according to the order used by <tt>MultiArrayView::operator[]</tt>.
1945  If the matrix is actually a vector, this is just the row or columns index.
1946  In case of a truely 2-dimensional matrix, the index can be converted to an
1947  index pair by calling <tt>MultiArrayView::scanOrderIndexToCoordinate()</tt>
1948 
1949  <b>Required Interface:</b>
1950 
1951  \code
1952  bool c = condition(a[0]);
1953  bool f = NumericTraits<T>::min() < a[0];
1954  \endcode
1955 
1956  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>><br>
1957  Namespace: vigra
1958  */
1959 template <class T, class C, class UnaryFunctor>
1960 int argMaxIf(MultiArrayView<2, T, C> const & a, UnaryFunctor condition)
1961 {
1962  T vopt = NumericTraits<T>::min();
1963  int best = -1;
1964  for(int k=0; k < a.size(); ++k)
1965  {
1966  if(condition(a[k]) && vopt < a[k])
1967  {
1968  vopt = a[k];
1969  best = k;
1970  }
1971  }
1972  return best;
1973 }
1974 
1975 /** Matrix point-wise power.
1976 */
1977 template <class T, class C>
1978 linalg::TemporaryMatrix<T> pow(MultiArrayView<2, T, C> const & v, T exponent)
1979 {
1980  linalg::TemporaryMatrix<T> t(v.shape());
1981  MultiArrayIndex m = rowCount(v), n = columnCount(v);
1982 
1983  for(MultiArrayIndex i = 0; i < n; ++i)
1984  for(MultiArrayIndex j = 0; j < m; ++j)
1985  t(j, i) = vigra::pow(v(j, i), exponent);
1986  return t;
1987 }
1988 
1989 template <class T>
1990 linalg::TemporaryMatrix<T> pow(linalg::TemporaryMatrix<T> const & v, T exponent)
1991 {
1992  linalg::TemporaryMatrix<T> & t = const_cast<linalg::TemporaryMatrix<T> &>(v);
1993  MultiArrayIndex m = rowCount(t), n = columnCount(t);
1994 
1995  for(MultiArrayIndex i = 0; i < n; ++i)
1996  for(MultiArrayIndex j = 0; j < m; ++j)
1997  t(j, i) = vigra::pow(t(j, i), exponent);
1998  return t;
1999 }
2000 
2001 template <class T, class C>
2002 linalg::TemporaryMatrix<T> pow(MultiArrayView<2, T, C> const & v, int exponent)
2003 {
2004  linalg::TemporaryMatrix<T> t(v.shape());
2005  MultiArrayIndex m = rowCount(v), n = columnCount(v);
2006 
2007  for(MultiArrayIndex i = 0; i < n; ++i)
2008  for(MultiArrayIndex j = 0; j < m; ++j)
2009  t(j, i) = vigra::pow(v(j, i), exponent);
2010  return t;
2011 }
2012 
2013 template <class T>
2014 linalg::TemporaryMatrix<T> pow(linalg::TemporaryMatrix<T> const & v, int exponent)
2015 {
2016  linalg::TemporaryMatrix<T> & t = const_cast<linalg::TemporaryMatrix<T> &>(v);
2017  MultiArrayIndex m = rowCount(t), n = columnCount(t);
2018 
2019  for(MultiArrayIndex i = 0; i < n; ++i)
2020  for(MultiArrayIndex j = 0; j < m; ++j)
2021  t(j, i) = vigra::pow(t(j, i), exponent);
2022  return t;
2023 }
2024 
2025 template <class C>
2026 linalg::TemporaryMatrix<int> pow(MultiArrayView<2, int, C> const & v, int exponent)
2027 {
2028  linalg::TemporaryMatrix<int> t(v.shape());
2029  MultiArrayIndex m = rowCount(v), n = columnCount(v);
2030 
2031  for(MultiArrayIndex i = 0; i < n; ++i)
2032  for(MultiArrayIndex j = 0; j < m; ++j)
2033  t(j, i) = (int)vigra::pow((double)v(j, i), exponent);
2034  return t;
2035 }
2036 
2037 inline
2038 linalg::TemporaryMatrix<int> pow(linalg::TemporaryMatrix<int> const & v, int exponent)
2039 {
2040  linalg::TemporaryMatrix<int> & t = const_cast<linalg::TemporaryMatrix<int> &>(v);
2041  MultiArrayIndex m = rowCount(t), n = columnCount(t);
2042 
2043  for(MultiArrayIndex i = 0; i < n; ++i)
2044  for(MultiArrayIndex j = 0; j < m; ++j)
2045  t(j, i) = (int)vigra::pow((double)t(j, i), exponent);
2046  return t;
2047 }
2048 
2049  /** Matrix point-wise sqrt. */
2050 template <class T, class C>
2051 linalg::TemporaryMatrix<T> sqrt(MultiArrayView<2, T, C> const & v);
2052  /** Matrix point-wise exp. */
2053 template <class T, class C>
2054 linalg::TemporaryMatrix<T> exp(MultiArrayView<2, T, C> const & v);
2055  /** Matrix point-wise log. */
2056 template <class T, class C>
2057 linalg::TemporaryMatrix<T> log(MultiArrayView<2, T, C> const & v);
2058  /** Matrix point-wise log10. */
2059 template <class T, class C>
2060 linalg::TemporaryMatrix<T> log10(MultiArrayView<2, T, C> const & v);
2061  /** Matrix point-wise sin. */
2062 template <class T, class C>
2063 linalg::TemporaryMatrix<T> sin(MultiArrayView<2, T, C> const & v);
2064  /** Matrix point-wise asin. */
2065 template <class T, class C>
2066 linalg::TemporaryMatrix<T> asin(MultiArrayView<2, T, C> const & v);
2067  /** Matrix point-wise cos. */
2068 template <class T, class C>
2069 linalg::TemporaryMatrix<T> cos(MultiArrayView<2, T, C> const & v);
2070  /** Matrix point-wise acos. */
2071 template <class T, class C>
2072 linalg::TemporaryMatrix<T> acos(MultiArrayView<2, T, C> const & v);
2073  /** Matrix point-wise tan. */
2074 template <class T, class C>
2075 linalg::TemporaryMatrix<T> tan(MultiArrayView<2, T, C> const & v);
2076  /** Matrix point-wise atan. */
2077 template <class T, class C>
2078 linalg::TemporaryMatrix<T> atan(MultiArrayView<2, T, C> const & v);
2079  /** Matrix point-wise round. */
2080 template <class T, class C>
2081 linalg::TemporaryMatrix<T> round(MultiArrayView<2, T, C> const & v);
2082  /** Matrix point-wise floor. */
2083 template <class T, class C>
2084 linalg::TemporaryMatrix<T> floor(MultiArrayView<2, T, C> const & v);
2085  /** Matrix point-wise ceil. */
2086 template <class T, class C>
2087 linalg::TemporaryMatrix<T> ceil(MultiArrayView<2, T, C> const & v);
2088  /** Matrix point-wise abs. */
2089 template <class T, class C>
2090 linalg::TemporaryMatrix<T> abs(MultiArrayView<2, T, C> const & v);
2091  /** Matrix point-wise square. */
2092 template <class T, class C>
2093 linalg::TemporaryMatrix<T> sq(MultiArrayView<2, T, C> const & v);
2094  /** Matrix point-wise sign. */
2095 template <class T, class C>
2096 linalg::TemporaryMatrix<T> sign(MultiArrayView<2, T, C> const & v);
2097 
2098 #define VIGRA_MATRIX_UNARY_FUNCTION(FUNCTION, NAMESPACE) \
2099 using NAMESPACE::FUNCTION; \
2100 template <class T, class C> \
2101 linalg::TemporaryMatrix<T> FUNCTION(MultiArrayView<2, T, C> const & v) \
2102 { \
2103  linalg::TemporaryMatrix<T> t(v.shape()); \
2104  MultiArrayIndex m = rowCount(v), n = columnCount(v); \
2105  \
2106  for(MultiArrayIndex i = 0; i < n; ++i) \
2107  for(MultiArrayIndex j = 0; j < m; ++j) \
2108  t(j, i) = NAMESPACE::FUNCTION(v(j, i)); \
2109  return t; \
2110 } \
2111  \
2112 template <class T> \
2113 linalg::TemporaryMatrix<T> FUNCTION(linalg::Matrix<T> const & v) \
2114 { \
2115  linalg::TemporaryMatrix<T> t(v.shape()); \
2116  MultiArrayIndex m = rowCount(v), n = columnCount(v); \
2117  \
2118  for(MultiArrayIndex i = 0; i < n; ++i) \
2119  for(MultiArrayIndex j = 0; j < m; ++j) \
2120  t(j, i) = NAMESPACE::FUNCTION(v(j, i)); \
2121  return t; \
2122 } \
2123  \
2124 template <class T> \
2125 linalg::TemporaryMatrix<T> FUNCTION(linalg::TemporaryMatrix<T> const & v) \
2126 { \
2127  linalg::TemporaryMatrix<T> & t = const_cast<linalg::TemporaryMatrix<T> &>(v); \
2128  MultiArrayIndex m = rowCount(t), n = columnCount(t); \
2129  \
2130  for(MultiArrayIndex i = 0; i < n; ++i) \
2131  for(MultiArrayIndex j = 0; j < m; ++j) \
2132  t(j, i) = NAMESPACE::FUNCTION(t(j, i)); \
2133  return v; \
2134 }\
2135 }\
2136 using linalg::FUNCTION;\
2137 namespace linalg {
2138 
2139 VIGRA_MATRIX_UNARY_FUNCTION(sqrt, std)
2140 VIGRA_MATRIX_UNARY_FUNCTION(exp, std)
2141 VIGRA_MATRIX_UNARY_FUNCTION(log, std)
2142 VIGRA_MATRIX_UNARY_FUNCTION(log10, std)
2143 VIGRA_MATRIX_UNARY_FUNCTION(sin, std)
2144 VIGRA_MATRIX_UNARY_FUNCTION(asin, std)
2145 VIGRA_MATRIX_UNARY_FUNCTION(cos, std)
2146 VIGRA_MATRIX_UNARY_FUNCTION(acos, std)
2147 VIGRA_MATRIX_UNARY_FUNCTION(tan, std)
2148 VIGRA_MATRIX_UNARY_FUNCTION(atan, std)
2149 VIGRA_MATRIX_UNARY_FUNCTION(round, vigra)
2150 VIGRA_MATRIX_UNARY_FUNCTION(floor, vigra)
2151 VIGRA_MATRIX_UNARY_FUNCTION(ceil, vigra)
2152 VIGRA_MATRIX_UNARY_FUNCTION(abs, vigra)
2153 VIGRA_MATRIX_UNARY_FUNCTION(sq, vigra)
2154 VIGRA_MATRIX_UNARY_FUNCTION(sign, vigra)
2155 
2156 #undef VIGRA_MATRIX_UNARY_FUNCTION
2157 
2158 //@}
2159 
2160 } // namespace linalg
2161 
2162 using linalg::RowMajor;
2163 using linalg::ColumnMajor;
2164 using linalg::Matrix;
2167 using linalg::transpose;
2168 using linalg::pointWise;
2169 using linalg::trace;
2170 using linalg::dot;
2171 using linalg::cross;
2172 using linalg::outer;
2173 using linalg::rowCount;
2174 using linalg::columnCount;
2175 using linalg::rowVector;
2176 using linalg::columnVector;
2177 using linalg::subVector;
2178 using linalg::isSymmetric;
2181 using linalg::argMin;
2182 using linalg::argMinIf;
2183 using linalg::argMax;
2184 using linalg::argMaxIf;
2185 
2186 /********************************************************/
2187 /* */
2188 /* NormTraits */
2189 /* */
2190 /********************************************************/
2191 
2192 template <class T, class ALLOC>
2193 struct NormTraits<Matrix<T, ALLOC> >
2194 : public NormTraits<MultiArray<2, T, ALLOC> >
2195 {
2196  typedef NormTraits<MultiArray<2, T, ALLOC> > BaseType;
2197  typedef Matrix<T, ALLOC> Type;
2198  typedef typename BaseType::SquaredNormType SquaredNormType;
2199  typedef typename BaseType::NormType NormType;
2200 };
2201 
2202 template <class T, class ALLOC>
2203 struct NormTraits<linalg::TemporaryMatrix<T, ALLOC> >
2204 : public NormTraits<Matrix<T, ALLOC> >
2205 {
2206  typedef NormTraits<Matrix<T, ALLOC> > BaseType;
2207  typedef linalg::TemporaryMatrix<T, ALLOC> Type;
2208  typedef typename BaseType::SquaredNormType SquaredNormType;
2209  typedef typename BaseType::NormType NormType;
2210 };
2211 
2212 } // namespace vigra
2213 
2214 namespace std {
2215 
2216 /** \addtogroup LinearAlgebraFunctions
2217  */
2218 //@{
2219 
2220  /** print a matrix \a m to the stream \a s.
2221 
2222  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
2223  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
2224  Namespace: std
2225  */
2226 template <class T, class C>
2227 ostream &
2228 operator<<(ostream & s, const vigra::MultiArrayView<2, T, C> &m)
2229 {
2232  ios::fmtflags flags = s.setf(ios::right | ios::fixed, ios::adjustfield | ios::floatfield);
2233  for(vigra::MultiArrayIndex j = 0; j < rows; ++j)
2234  {
2235  for(vigra::MultiArrayIndex i = 0; i < cols; ++i)
2236  {
2237  s << m(j, i) << " ";
2238  }
2239  s << endl;
2240  }
2241  s.setf(flags);
2242  return s;
2243 }
2244 
2245 //@}
2246 
2247 } // namespace std
2248 
2249 namespace vigra {
2250 
2251 namespace linalg {
2252 
2253 namespace detail {
2254 
2255 template <class T1, class C1, class T2, class C2, class T3, class C3>
2256 void
2257 columnStatisticsImpl(MultiArrayView<2, T1, C1> const & A,
2258  MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & sumOfSquaredDifferences)
2259 {
2260  MultiArrayIndex m = rowCount(A);
2262  vigra_precondition(1 == rowCount(mean) && n == columnCount(mean) &&
2263  1 == rowCount(sumOfSquaredDifferences) && n == columnCount(sumOfSquaredDifferences),
2264  "columnStatistics(): Shape mismatch between input and output.");
2265 
2266  // West's algorithm for incremental variance computation
2267  mean.init(NumericTraits<T2>::zero());
2268  sumOfSquaredDifferences.init(NumericTraits<T3>::zero());
2269 
2270  for(MultiArrayIndex k=0; k<m; ++k)
2271  {
2272  Matrix<T2> t = rowVector(A, k) - mean;
2273  typename NumericTraits<T2>::RealPromote f = 1.0 / (k + 1.0),
2274  f1 = 1.0 - f;
2275  mean += f*t;
2276  sumOfSquaredDifferences += f1*sq(t);
2277  }
2278 }
2279 
2280 template <class T1, class C1, class T2, class C2, class T3, class C3>
2281 void
2282 columnStatistics2PassImpl(MultiArrayView<2, T1, C1> const & A,
2283  MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & sumOfSquaredDifferences)
2284 {
2285  MultiArrayIndex m = rowCount(A);
2287  vigra_precondition(1 == rowCount(mean) && n == columnCount(mean) &&
2288  1 == rowCount(sumOfSquaredDifferences) && n == columnCount(sumOfSquaredDifferences),
2289  "columnStatistics(): Shape mismatch between input and output.");
2290 
2291  // two-pass algorithm for incremental variance computation
2292  mean.init(NumericTraits<T2>::zero());
2293  for(MultiArrayIndex k=0; k<m; ++k)
2294  {
2295  mean += rowVector(A, k);
2296  }
2297  mean /= (double)m;
2298 
2299  sumOfSquaredDifferences.init(NumericTraits<T3>::zero());
2300  for(MultiArrayIndex k=0; k<m; ++k)
2301  {
2302  sumOfSquaredDifferences += sq(rowVector(A, k) - mean);
2303  }
2304 }
2305 
2306 } // namespace detail
2307 
2308 /** \addtogroup LinearAlgebraFunctions
2309  */
2310 //@{
2311  /** Compute statistics of every column of matrix \a A.
2312 
2313  The result matrices must be row vectors with as many columns as \a A.
2314 
2315  <b> Declarations:</b>
2316 
2317  compute only the mean:
2318  \code
2319  namespace vigra { namespace linalg {
2320  template <class T1, class C1, class T2, class C2>
2321  void
2322  columnStatistics(MultiArrayView<2, T1, C1> const & A,
2323  MultiArrayView<2, T2, C2> & mean);
2324  } }
2325  \endcode
2326 
2327  compute mean and standard deviation:
2328  \code
2329  namespace vigra { namespace linalg {
2330  template <class T1, class C1, class T2, class C2, class T3, class C3>
2331  void
2332  columnStatistics(MultiArrayView<2, T1, C1> const & A,
2333  MultiArrayView<2, T2, C2> & mean,
2334  MultiArrayView<2, T3, C3> & stdDev);
2335  } }
2336  \endcode
2337 
2338  compute mean, standard deviation, and norm:
2339  \code
2340  namespace vigra { namespace linalg {
2341  template <class T1, class C1, class T2, class C2, class T3, class C3, class T4, class C4>
2342  void
2343  columnStatistics(MultiArrayView<2, T1, C1> const & A,
2344  MultiArrayView<2, T2, C2> & mean,
2345  MultiArrayView<2, T3, C3> & stdDev,
2346  MultiArrayView<2, T4, C4> & norm);
2347  } }
2348  \endcode
2349 
2350  <b> Usage:</b>
2351 
2352  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
2353  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
2354  Namespaces: vigra and vigra::linalg
2355 
2356  \code
2357  Matrix A(rows, columns);
2358  .. // fill A
2359  Matrix columnMean(1, columns), columnStdDev(1, columns), columnNorm(1, columns);
2360 
2361  columnStatistics(A, columnMean, columnStdDev, columnNorm);
2362 
2363  \endcode
2364  */
2365 doxygen_overloaded_function(template <...> void columnStatistics)
2366 
2367 template <class T1, class C1, class T2, class C2>
2368 void
2369 columnStatistics(MultiArrayView<2, T1, C1> const & A,
2370  MultiArrayView<2, T2, C2> & mean)
2371 {
2372  MultiArrayIndex m = rowCount(A);
2374  vigra_precondition(1 == rowCount(mean) && n == columnCount(mean),
2375  "columnStatistics(): Shape mismatch between input and output.");
2376 
2377  mean.init(NumericTraits<T2>::zero());
2378 
2379  for(MultiArrayIndex k=0; k<m; ++k)
2380  {
2381  mean += rowVector(A, k);
2382  }
2383  mean /= T2(m);
2384 }
2385 
2386 template <class T1, class C1, class T2, class C2, class T3, class C3>
2387 void
2388 columnStatistics(MultiArrayView<2, T1, C1> const & A,
2389  MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & stdDev)
2390 {
2391  detail::columnStatisticsImpl(A, mean, stdDev);
2392 
2393  if(rowCount(A) > 1)
2394  stdDev = sqrt(stdDev / T3(rowCount(A) - 1.0));
2395 }
2396 
2397 template <class T1, class C1, class T2, class C2, class T3, class C3, class T4, class C4>
2398 void
2399 columnStatistics(MultiArrayView<2, T1, C1> const & A,
2400  MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & stdDev, MultiArrayView<2, T4, C4> & norm)
2401 {
2402  MultiArrayIndex m = rowCount(A);
2404  vigra_precondition(1 == rowCount(mean) && n == columnCount(mean) &&
2405  1 == rowCount(stdDev) && n == columnCount(stdDev) &&
2406  1 == rowCount(norm) && n == columnCount(norm),
2407  "columnStatistics(): Shape mismatch between input and output.");
2408 
2409  detail::columnStatisticsImpl(A, mean, stdDev);
2410  norm = sqrt(stdDev + T2(m) * sq(mean));
2411  stdDev = sqrt(stdDev / T3(m - 1.0));
2412 }
2413 
2414  /** Compute statistics of every row of matrix \a A.
2415 
2416  The result matrices must be column vectors with as many rows as \a A.
2417 
2418  <b> Declarations:</b>
2419 
2420  compute only the mean:
2421  \code
2422  namespace vigra { namespace linalg {
2423  template <class T1, class C1, class T2, class C2>
2424  void
2425  rowStatistics(MultiArrayView<2, T1, C1> const & A,
2426  MultiArrayView<2, T2, C2> & mean);
2427  } }
2428  \endcode
2429 
2430  compute mean and standard deviation:
2431  \code
2432  namespace vigra { namespace linalg {
2433  template <class T1, class C1, class T2, class C2, class T3, class C3>
2434  void
2435  rowStatistics(MultiArrayView<2, T1, C1> const & A,
2436  MultiArrayView<2, T2, C2> & mean,
2437  MultiArrayView<2, T3, C3> & stdDev);
2438  } }
2439  \endcode
2440 
2441  compute mean, standard deviation, and norm:
2442  \code
2443  namespace vigra { namespace linalg {
2444  template <class T1, class C1, class T2, class C2, class T3, class C3, class T4, class C4>
2445  void
2446  rowStatistics(MultiArrayView<2, T1, C1> const & A,
2447  MultiArrayView<2, T2, C2> & mean,
2448  MultiArrayView<2, T3, C3> & stdDev,
2449  MultiArrayView<2, T4, C4> & norm);
2450  } }
2451  \endcode
2452 
2453  <b> Usage:</b>
2454 
2455  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
2456  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
2457  Namespaces: vigra and vigra::linalg
2458 
2459  \code
2460  Matrix A(rows, columns);
2461  .. // fill A
2462  Matrix rowMean(rows, 1), rowStdDev(rows, 1), rowNorm(rows, 1);
2463 
2464  rowStatistics(a, rowMean, rowStdDev, rowNorm);
2465 
2466  \endcode
2467  */
2468 doxygen_overloaded_function(template <...> void rowStatistics)
2469 
2470 template <class T1, class C1, class T2, class C2>
2471 void
2472 rowStatistics(MultiArrayView<2, T1, C1> const & A,
2473  MultiArrayView<2, T2, C2> & mean)
2474 {
2475  vigra_precondition(1 == columnCount(mean) && rowCount(A) == rowCount(mean),
2476  "rowStatistics(): Shape mismatch between input and output.");
2477  MultiArrayView<2, T2, StridedArrayTag> tm = transpose(mean);
2478  columnStatistics(transpose(A), tm);
2479 }
2480 
2481 template <class T1, class C1, class T2, class C2, class T3, class C3>
2482 void
2483 rowStatistics(MultiArrayView<2, T1, C1> const & A,
2484  MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & stdDev)
2485 {
2486  vigra_precondition(1 == columnCount(mean) && rowCount(A) == rowCount(mean) &&
2487  1 == columnCount(stdDev) && rowCount(A) == rowCount(stdDev),
2488  "rowStatistics(): Shape mismatch between input and output.");
2489  MultiArrayView<2, T2, StridedArrayTag> tm = transpose(mean);
2490  MultiArrayView<2, T3, StridedArrayTag> ts = transpose(stdDev);
2491  columnStatistics(transpose(A), tm, ts);
2492 }
2493 
2494 template <class T1, class C1, class T2, class C2, class T3, class C3, class T4, class C4>
2495 void
2496 rowStatistics(MultiArrayView<2, T1, C1> const & A,
2497  MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & stdDev, MultiArrayView<2, T4, C4> & norm)
2498 {
2499  vigra_precondition(1 == columnCount(mean) && rowCount(A) == rowCount(mean) &&
2500  1 == columnCount(stdDev) && rowCount(A) == rowCount(stdDev) &&
2501  1 == columnCount(norm) && rowCount(A) == rowCount(norm),
2502  "rowStatistics(): Shape mismatch between input and output.");
2503  MultiArrayView<2, T2, StridedArrayTag> tm = transpose(mean);
2504  MultiArrayView<2, T3, StridedArrayTag> ts = transpose(stdDev);
2505  MultiArrayView<2, T4, StridedArrayTag> tn = transpose(norm);
2506  columnStatistics(transpose(A), tm, ts, tn);
2507 }
2508 
2509 namespace detail {
2510 
2511 template <class T1, class C1, class U, class T2, class C2, class T3, class C3>
2512 void updateCovarianceMatrix(MultiArrayView<2, T1, C1> const & features,
2513  U & count, MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & covariance)
2514 {
2515  MultiArrayIndex n = std::max(rowCount(features), columnCount(features));
2516  vigra_precondition(std::min(rowCount(features), columnCount(features)) == 1,
2517  "updateCovarianceMatrix(): Features must be a row or column vector.");
2518  vigra_precondition(mean.shape() == features.shape(),
2519  "updateCovarianceMatrix(): Shape mismatch between feature vector and mean vector.");
2520  vigra_precondition(n == rowCount(covariance) && n == columnCount(covariance),
2521  "updateCovarianceMatrix(): Shape mismatch between feature vector and covariance matrix.");
2522 
2523  // West's algorithm for incremental covariance matrix computation
2524  Matrix<T2> t = features - mean;
2525  ++count;
2526  double f = 1.0 / count,
2527  f1 = 1.0 - f;
2528  mean += f*t;
2529 
2530  if(rowCount(features) == 1) // update column covariance from current row
2531  {
2532  for(MultiArrayIndex k=0; k<n; ++k)
2533  {
2534  covariance(k, k) += f1*sq(t(0, k));
2535  for(MultiArrayIndex l=k+1; l<n; ++l)
2536  {
2537  covariance(k, l) += f1*t(0, k)*t(0, l);
2538  covariance(l, k) = covariance(k, l);
2539  }
2540  }
2541  }
2542  else // update row covariance from current column
2543  {
2544  for(MultiArrayIndex k=0; k<n; ++k)
2545  {
2546  covariance(k, k) += f1*sq(t(k, 0));
2547  for(MultiArrayIndex l=k+1; l<n; ++l)
2548  {
2549  covariance(k, l) += f1*t(k, 0)*t(l, 0);
2550  covariance(l, k) = covariance(k, l);
2551  }
2552  }
2553  }
2554 }
2555 
2556 } // namespace detail
2557 
2558  /*! Compute the covariance matrix between the columns of a matrix \a features.
2559 
2560  The result matrix \a covariance must by a square matrix with as many rows and
2561  columns as the number of columns in matrix \a features.
2562 
2563  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>><br>
2564  Namespace: vigra
2565  */
2566 template <class T1, class C1, class T2, class C2>
2568  MultiArrayView<2, T2, C2> & covariance)
2569 {
2570  MultiArrayIndex m = rowCount(features), n = columnCount(features);
2571  vigra_precondition(n == rowCount(covariance) && n == columnCount(covariance),
2572  "covarianceMatrixOfColumns(): Shape mismatch between feature matrix and covariance matrix.");
2573  MultiArrayIndex count = 0;
2574  Matrix<T2> means(1, n);
2575  covariance.init(NumericTraits<T2>::zero());
2576  for(MultiArrayIndex k=0; k<m; ++k)
2577  detail::updateCovarianceMatrix(rowVector(features, k), count, means, covariance);
2578  covariance /= T2(m - 1);
2579 }
2580 
2581  /*! Compute the covariance matrix between the columns of a matrix \a features.
2582 
2583  The result is returned as a square temporary matrix with as many rows and
2584  columns as the number of columns in matrix \a features.
2585 
2586  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>><br>
2587  Namespace: vigra
2588  */
2589 template <class T, class C>
2590 TemporaryMatrix<T>
2592 {
2593  TemporaryMatrix<T> res(columnCount(features), columnCount(features));
2594  covarianceMatrixOfColumns(features, res);
2595  return res;
2596 }
2597 
2598  /*! Compute the covariance matrix between the rows of a matrix \a features.
2599 
2600  The result matrix \a covariance must by a square matrix with as many rows and
2601  columns as the number of rows in matrix \a features.
2602 
2603  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>><br>
2604  Namespace: vigra
2605  */
2606 template <class T1, class C1, class T2, class C2>
2608  MultiArrayView<2, T2, C2> & covariance)
2609 {
2610  MultiArrayIndex m = rowCount(features), n = columnCount(features);
2611  vigra_precondition(m == rowCount(covariance) && m == columnCount(covariance),
2612  "covarianceMatrixOfRows(): Shape mismatch between feature matrix and covariance matrix.");
2613  MultiArrayIndex count = 0;
2614  Matrix<T2> means(m, 1);
2615  covariance.init(NumericTraits<T2>::zero());
2616  for(MultiArrayIndex k=0; k<n; ++k)
2617  detail::updateCovarianceMatrix(columnVector(features, k), count, means, covariance);
2618  covariance /= T2(m - 1);
2619 }
2620 
2621  /*! Compute the covariance matrix between the rows of a matrix \a features.
2622 
2623  The result is returned as a square temporary matrix with as many rows and
2624  columns as the number of rows in matrix \a features.
2625 
2626  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>><br>
2627  Namespace: vigra
2628  */
2629 template <class T, class C>
2630 TemporaryMatrix<T>
2632 {
2633  TemporaryMatrix<T> res(rowCount(features), rowCount(features));
2634  covarianceMatrixOfRows(features, res);
2635  return res;
2636 }
2637 
2638 enum DataPreparationGoals { ZeroMean = 1, UnitVariance = 2, UnitNorm = 4 };
2639 
2640 inline DataPreparationGoals operator|(DataPreparationGoals l, DataPreparationGoals r)
2641 {
2642  return DataPreparationGoals(int(l) | int(r));
2643 }
2644 
2645 namespace detail {
2646 
2647 template <class T, class C1, class C2, class C3, class C4>
2648 void
2649 prepareDataImpl(const MultiArrayView<2, T, C1> & A,
2650  MultiArrayView<2, T, C2> & res, MultiArrayView<2, T, C3> & offset, MultiArrayView<2, T, C4> & scaling,
2651  DataPreparationGoals goals)
2652 {
2653  MultiArrayIndex m = rowCount(A);
2655  vigra_precondition(A.shape() == res.shape() &&
2656  n == columnCount(offset) && 1 == rowCount(offset) &&
2657  n == columnCount(scaling) && 1 == rowCount(scaling),
2658  "prepareDataImpl(): Shape mismatch between input and output.");
2659 
2660  if(!goals)
2661  {
2662  res = A;
2663  offset.init(NumericTraits<T>::zero());
2664  scaling.init(NumericTraits<T>::one());
2665  return;
2666  }
2667 
2668  bool zeroMean = (goals & ZeroMean) != 0;
2669  bool unitVariance = (goals & UnitVariance) != 0;
2670  bool unitNorm = (goals & UnitNorm) != 0;
2671 
2672  vigra_precondition(!(unitVariance && unitNorm),
2673  "prepareDataImpl(): Unit variance and unit norm cannot be achieved at the same time.");
2674 
2675  Matrix<T> mean(1, n), sumOfSquaredDifferences(1, n);
2676  detail::columnStatisticsImpl(A, mean, sumOfSquaredDifferences);
2677 
2678  for(MultiArrayIndex k=0; k<n; ++k)
2679  {
2680  T stdDev = std::sqrt(sumOfSquaredDifferences(0, k) / T(m-1));
2681  if(closeAtTolerance(stdDev / mean(0,k), NumericTraits<T>::zero()))
2682  stdDev = NumericTraits<T>::zero();
2683  if(zeroMean && stdDev > NumericTraits<T>::zero())
2684  {
2685  columnVector(res, k) = columnVector(A, k) - mean(0,k);
2686  offset(0, k) = mean(0, k);
2687  mean(0, k) = NumericTraits<T>::zero();
2688  }
2689  else
2690  {
2691  columnVector(res, k) = columnVector(A, k);
2692  offset(0, k) = NumericTraits<T>::zero();
2693  }
2694 
2695  T norm = mean(0,k) == NumericTraits<T>::zero()
2696  ? std::sqrt(sumOfSquaredDifferences(0, k))
2697  : std::sqrt(sumOfSquaredDifferences(0, k) + T(m) * sq(mean(0,k)));
2698  if(unitNorm && norm > NumericTraits<T>::zero())
2699  {
2700  columnVector(res, k) /= norm;
2701  scaling(0, k) = NumericTraits<T>::one() / norm;
2702  }
2703  else if(unitVariance && stdDev > NumericTraits<T>::zero())
2704  {
2705  columnVector(res, k) /= stdDev;
2706  scaling(0, k) = NumericTraits<T>::one() / stdDev;
2707  }
2708  else
2709  {
2710  scaling(0, k) = NumericTraits<T>::one();
2711  }
2712  }
2713 }
2714 
2715 } // namespace detail
2716 
2717  /*! Standardize the columns of a matrix according to given <tt>DataPreparationGoals</tt>.
2718 
2719  For every column of the matrix \a A, this function computes mean,
2720  standard deviation, and norm. It then applies a linear transformation to the values of
2721  the column according to these statistics and the given <tt>DataPreparationGoals</tt>.
2722  The result is returned in matrix \a res which must have the same size as \a A.
2723  Optionally, the transformation applied can also be returned in the matrices \a offset
2724  and \a scaling (see below for an example how these matrices can be used to standardize
2725  more data according to the same transformation).
2726 
2727  The following <tt>DataPreparationGoals</tt> are supported:
2728 
2729  <DL>
2730  <DT><tt>ZeroMean</tt><DD> Subtract the column mean form every column if the values in the column are not constant.
2731  Do nothing in a constant column.
2732  <DT><tt>UnitVariance</tt><DD> Divide by the column standard deviation if the values in the column are not constant.
2733  Do nothing in a constant column.
2734  <DT><tt>UnitNorm</tt><DD> Divide by the column norm if it is non-zero.
2735  <DT><tt>ZeroMean | UnitVariance</tt><DD> First subtact the mean and then divide by the standard deviation, unless the
2736  column is constant (in which case the column remains unchanged).
2737  <DT><tt>ZeroMean | UnitNorm</tt><DD> If the column is non-constant, subtract the mean. Then divide by the norm
2738  of the result if the norm is non-zero.
2739  </DL>
2740 
2741  <b> Declarations:</b>
2742 
2743  Standardize the matrix and return the parameters of the linear transformation.
2744  The matrices \a offset and \a scaling must be row vectors with as many columns as \a A.
2745  \code
2746  namespace vigra { namespace linalg {
2747  template <class T, class C1, class C2, class C3, class C4>
2748  void
2749  prepareColumns(MultiArrayView<2, T, C1> const & A,
2750  MultiArrayView<2, T, C2> & res,
2751  MultiArrayView<2, T, C3> & offset,
2752  MultiArrayView<2, T, C4> & scaling,
2753  DataPreparationGoals goals = ZeroMean | UnitVariance);
2754  } }
2755  \endcode
2756 
2757  Only standardize the matrix.
2758  \code
2759  namespace vigra { namespace linalg {
2760  template <class T, class C1, class C2>
2761  void
2762  prepareColumns(MultiArrayView<2, T, C1> const & A,
2763  MultiArrayView<2, T, C2> & res,
2764  DataPreparationGoals goals = ZeroMean | UnitVariance);
2765  } }
2766  \endcode
2767 
2768  <b> Usage:</b>
2769 
2770  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
2771  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
2772  Namespaces: vigra and vigra::linalg
2773 
2774  \code
2775  Matrix A(rows, columns);
2776  .. // fill A
2777  Matrix standardizedA(rows, columns), offset(1, columns), scaling(1, columns);
2778 
2779  prepareColumns(A, standardizedA, offset, scaling, ZeroMean | UnitNorm);
2780 
2781  // use offset and scaling to prepare additional data according to the same transformation
2782  Matrix newData(nrows, columns);
2783 
2784  Matrix standardizedNewData = (newData - repeatMatrix(offset, nrows, 1)) * pointWise(repeatMatrix(scaling, nrows, 1));
2785 
2786  \endcode
2787  */
2788 doxygen_overloaded_function(template <...> void prepareColumns)
2789 
2790 template <class T, class C1, class C2, class C3, class C4>
2791 inline void
2792 prepareColumns(MultiArrayView<2, T, C1> const & A,
2793  MultiArrayView<2, T, C2> & res, MultiArrayView<2, T, C3> & offset, MultiArrayView<2, T, C4> & scaling,
2794  DataPreparationGoals goals = ZeroMean | UnitVariance)
2795 {
2796  detail::prepareDataImpl(A, res, offset, scaling, goals);
2797 }
2798 
2799 template <class T, class C1, class C2>
2800 inline void
2801 prepareColumns(MultiArrayView<2, T, C1> const & A, MultiArrayView<2, T, C2> & res,
2802  DataPreparationGoals goals = ZeroMean | UnitVariance)
2803 {
2804  Matrix<T> offset(1, columnCount(A)), scaling(1, columnCount(A));
2805  detail::prepareDataImpl(A, res, offset, scaling, goals);
2806 }
2807 
2808  /*! Standardize the rows of a matrix according to given <tt>DataPreparationGoals</tt>.
2809 
2810  This algorithm works in the same way as \ref prepareColumns() (see there for detailed
2811  documentation), but is applied to the rows of the matrix \a A instead. Accordingly, the
2812  matrices holding the parameters of the linear transformation must be column vectors
2813  with as many rows as \a A.
2814 
2815  <b> Declarations:</b>
2816 
2817  Standardize the matrix and return the parameters of the linear transformation.
2818  The matrices \a offset and \a scaling must be column vectors
2819  with as many rows as \a A.
2820 
2821  \code
2822  namespace vigra { namespace linalg {
2823  template <class T, class C1, class C2, class C3, class C4>
2824  void
2825  prepareRows(MultiArrayView<2, T, C1> const & A,
2826  MultiArrayView<2, T, C2> & res,
2827  MultiArrayView<2, T, C3> & offset,
2828  MultiArrayView<2, T, C4> & scaling,
2829  DataPreparationGoals goals = ZeroMean | UnitVariance)´;
2830  } }
2831  \endcode
2832 
2833  Only standardize the matrix.
2834  \code
2835  namespace vigra { namespace linalg {
2836  template <class T, class C1, class C2>
2837  void
2838  prepareRows(MultiArrayView<2, T, C1> const & A,
2839  MultiArrayView<2, T, C2> & res,
2840  DataPreparationGoals goals = ZeroMean | UnitVariance);
2841  } }
2842  \endcode
2843 
2844  <b> Usage:</b>
2845 
2846  <b>\#include</b> <<a href="matrix_8hxx-source.html">vigra/matrix.hxx</a>> or<br>
2847  <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
2848  Namespaces: vigra and vigra::linalg
2849 
2850  \code
2851  Matrix A(rows, columns);
2852  .. // fill A
2853  Matrix standardizedA(rows, columns), offset(rows, 1), scaling(rows, 1);
2854 
2855  prepareRows(A, standardizedA, offset, scaling, ZeroMean | UnitNorm);
2856 
2857  // use offset and scaling to prepare additional data according to the same transformation
2858  Matrix newData(rows, ncolumns);
2859 
2860  Matrix standardizedNewData = (newData - repeatMatrix(offset, 1, ncolumns)) * pointWise(repeatMatrix(scaling, 1, ncolumns));
2861 
2862  \endcode
2863 */
2864 doxygen_overloaded_function(template <...> void prepareRows)
2865 
2866 template <class T, class C1, class C2, class C3, class C4>
2867 inline void
2868 prepareRows(MultiArrayView<2, T, C1> const & A,
2869  MultiArrayView<2, T, C2> & res, MultiArrayView<2, T, C3> & offset, MultiArrayView<2, T, C4> & scaling,
2870  DataPreparationGoals goals = ZeroMean | UnitVariance)
2871 {
2872  MultiArrayView<2, T, StridedArrayTag> tr = transpose(res), to = transpose(offset), ts = transpose(scaling);
2873  detail::prepareDataImpl(transpose(A), tr, to, ts, goals);
2874 }
2875 
2876 template <class T, class C1, class C2>
2877 inline void
2878 prepareRows(MultiArrayView<2, T, C1> const & A, MultiArrayView<2, T, C2> & res,
2879  DataPreparationGoals goals = ZeroMean | UnitVariance)
2880 {
2881  MultiArrayView<2, T, StridedArrayTag> tr = transpose(res);
2882  Matrix<T> offset(rowCount(A), 1), scaling(rowCount(A), 1);
2883  detail::prepareDataImpl(transpose(A), tr, offset, scaling, goals);
2884 }
2885 
2886 //@}
2887 
2888 } // namespace linalg
2889 
2892 using linalg::rowStatistics;
2893 using linalg::prepareRows;
2894 using linalg::ZeroMean;
2895 using linalg::UnitVariance;
2896 using linalg::UnitNorm;
2897 
2898 } // namespace vigra
2899 
2900 
2901 
2902 #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.7.1 (Thu Jun 14 2012)