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

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

html generated using doxygen and Python
vigra 1.9.0 (Tue Oct 22 2013)