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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.8.0 (Wed Sep 26 2012)