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

matrix.hxx VIGRA

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, StridedArrayTag> 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  /** calculate the outer product of a TinyVector with itself.
1526  The result is returned as a temporary matrix.
1527 
1528  <b>\#include</b> <vigra/matrix.hxx> or<br>
1529  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1530  Namespaces: vigra and vigra::linalg
1531  */
1532 template <class T, int N>
1533 TemporaryMatrix<T>
1535 {
1536  TemporaryMatrix<T> ret(N, N);
1537 
1538  for(MultiArrayIndex i = 0; i < N; ++i)
1539  for(MultiArrayIndex j = 0; j < N; ++j)
1540  ret(j, i) = x[j] * x[i];
1541  return ret;
1542 }
1543 
1544 template <class T>
1545 class PointWise
1546 {
1547  public:
1548  T const & t;
1549 
1550  PointWise(T const & it)
1551  : t(it)
1552  {}
1553 };
1554 
1555 template <class T>
1556 PointWise<T> pointWise(T const & t)
1557 {
1558  return PointWise<T>(t);
1559 }
1560 
1561 
1562  /** multiply matrix \a a with scalar \a b.
1563  The result is written into \a r. \a a and \a r must have the same shape.
1564 
1565  <b>\#include</b> <vigra/matrix.hxx> or<br>
1566  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1567  Namespace: vigra::linalg
1568  */
1569 template <class T, class C1, class C2>
1571 {
1572  const MultiArrayIndex rows = rowCount(a);
1573  const MultiArrayIndex cols = columnCount(a);
1574  vigra_precondition(rows == rowCount(r) && cols == columnCount(r),
1575  "smul(): Matrix sizes must agree.");
1576 
1577  for(MultiArrayIndex i = 0; i < cols; ++i)
1578  for(MultiArrayIndex j = 0; j < rows; ++j)
1579  r(j, i) = a(j, i) * b;
1580 }
1581 
1582  /** multiply scalar \a a with matrix \a b.
1583  The result is written into \a r. \a b and \a r must have the same shape.
1584 
1585  <b>\#include</b> <vigra/matrix.hxx> or<br>
1586  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1587  Namespace: vigra::linalg
1588  */
1589 template <class T, class C2, class C3>
1591 {
1592  smul(b, a, r);
1593 }
1594 
1595  /** perform matrix multiplication of matrices \a a and \a b.
1596  The result is written into \a r. The three matrices must have matching shapes.
1597 
1598  <b>\#include</b> <vigra/matrix.hxx> or<br>
1599  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1600  Namespace: vigra::linalg
1601  */
1602 template <class T, class C1, class C2, class C3>
1605 {
1606  const MultiArrayIndex rrows = rowCount(r);
1607  const MultiArrayIndex rcols = columnCount(r);
1608  const MultiArrayIndex acols = columnCount(a);
1609  vigra_precondition(rrows == rowCount(a) && rcols == columnCount(b) && acols == rowCount(b),
1610  "mmul(): Matrix shapes must agree.");
1611 
1612  // order of loops ensures that inner loop goes down columns
1613  for(MultiArrayIndex i = 0; i < rcols; ++i)
1614  {
1615  for(MultiArrayIndex j = 0; j < rrows; ++j)
1616  r(j, i) = a(j, 0) * b(0, i);
1617  for(MultiArrayIndex k = 1; k < acols; ++k)
1618  for(MultiArrayIndex j = 0; j < rrows; ++j)
1619  r(j, i) += a(j, k) * b(k, i);
1620  }
1621 }
1622 
1623  /** perform matrix multiplication of matrices \a a and \a b.
1624  \a a and \a b must have matching shapes.
1625  The result is returned as a temporary matrix.
1626 
1627  <b>\#include</b> <vigra/matrix.hxx> or<br>
1628  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1629  Namespace: vigra::linalg
1630  */
1631 template <class T, class C1, class C2>
1632 inline TemporaryMatrix<T>
1634 {
1635  TemporaryMatrix<T> ret(rowCount(a), columnCount(b));
1636  mmul(a, b, ret);
1637  return ret;
1638 }
1639 
1640  /** multiply two matrices \a a and \a b pointwise.
1641  The result is written into \a r. All three matrices must have the same shape.
1642 
1643  <b>\#include</b> <vigra/matrix.hxx> or<br>
1644  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1645  Namespace: vigra::linalg
1646  */
1647 template <class T, class C1, class C2, class C3>
1650 {
1651  const MultiArrayIndex rrows = rowCount(r);
1652  const MultiArrayIndex rcols = columnCount(r);
1653  vigra_precondition(rrows == rowCount(a) && rcols == columnCount(a) &&
1654  rrows == rowCount(b) && rcols == columnCount(b),
1655  "pmul(): Matrix shapes must agree.");
1656 
1657  for(MultiArrayIndex i = 0; i < rcols; ++i) {
1658  for(MultiArrayIndex j = 0; j < rrows; ++j) {
1659  r(j, i) = a(j, i) * b(j, i);
1660  }
1661  }
1662 }
1663 
1664  /** multiply matrices \a a and \a b pointwise.
1665  \a a and \a b must have matching shapes.
1666  The result is returned as a temporary matrix.
1667 
1668  <b>\#include</b> <vigra/matrix.hxx> or<br>
1669  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1670  Namespace: vigra::linalg
1671  */
1672 template <class T, class C1, class C2>
1673 inline TemporaryMatrix<T>
1675 {
1676  TemporaryMatrix<T> ret(a.shape());
1677  pmul(a, b, ret);
1678  return ret;
1679 }
1680 
1681  /** multiply matrices \a a and \a b pointwise.
1682  \a a and \a b must have matching shapes.
1683  The result is returned as a temporary matrix.
1684 
1685  Usage:
1686 
1687  \code
1688  Matrix<double> a(m,n), b(m,n);
1689 
1690  Matrix<double> c = a * pointWise(b);
1691  // is equivalent to
1692  // Matrix<double> c = pmul(a, b);
1693  \endcode
1694 
1695  <b>\#include</b> <vigra/matrix.hxx> or<br>
1696  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1697  Namespace: vigra::linalg
1698  */
1699 template <class T, class C, class U>
1700 inline TemporaryMatrix<T>
1701 operator*(const MultiArrayView<2, T, C> &a, PointWise<U> b)
1702 {
1703  return pmul(a, b.t);
1704 }
1705 
1706  /** multiply matrix \a a with scalar \a b.
1707  The result is returned as a temporary matrix.
1708 
1709  <b>\#include</b> <vigra/matrix.hxx> or<br>
1710  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1711  Namespace: vigra::linalg
1712  */
1713 template <class T, class C>
1714 inline TemporaryMatrix<T>
1716 {
1717  return TemporaryMatrix<T>(a) *= b;
1718 }
1719 
1720 template <class T>
1721 inline TemporaryMatrix<T>
1722 operator*(const TemporaryMatrix<T> &a, T b)
1723 {
1724  return const_cast<TemporaryMatrix<T> &>(a) *= b;
1725 }
1726 
1727  /** multiply scalar \a a with matrix \a b.
1728  The result is returned as a temporary matrix.
1729 
1730  <b>\#include</b> <vigra/matrix.hxx> or<br>
1731  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1732  Namespace: vigra::linalg
1733  */
1734 template <class T, class C>
1735 inline TemporaryMatrix<T>
1737 {
1738  return TemporaryMatrix<T>(b) *= a;
1739 }
1740 
1741 template <class T>
1742 inline TemporaryMatrix<T>
1743 operator*(T a, const TemporaryMatrix<T> &b)
1744 {
1745  return const_cast<TemporaryMatrix<T> &>(b) *= a;
1746 }
1747 
1748  /** multiply matrix \a a with TinyVector \a b.
1749  \a a must be of size <tt>N x N</tt>. Vector \a b and the result
1750  vector are interpreted as column vectors.
1751 
1752  <b>\#include</b> <vigra/matrix.hxx> or<br>
1753  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1754  Namespace: vigra::linalg
1755  */
1756 template <class T, class A, int N, class DATA, class DERIVED>
1757 TinyVector<T, N>
1759 {
1760  vigra_precondition(N == rowCount(a) && N == columnCount(a),
1761  "operator*(Matrix, TinyVector): Shape mismatch.");
1762 
1763  TinyVector<T, N> res = TinyVectorView<T, N>(&a(0,0)) * b[0];
1764  for(MultiArrayIndex i = 1; i < N; ++i)
1765  res += TinyVectorView<T, N>(&a(0,i)) * b[i];
1766  return res;
1767 }
1768 
1769  /** multiply TinyVector \a a with matrix \a b.
1770  \a b must be of size <tt>N x N</tt>. Vector \a a and the result
1771  vector are interpreted as row vectors.
1772 
1773  <b>\#include</b> <vigra/matrix.hxx> or<br>
1774  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1775  Namespace: vigra::linalg
1776  */
1777 template <class T, int N, class DATA, class DERIVED, class A>
1780 {
1781  vigra_precondition(N == rowCount(b) && N == columnCount(b),
1782  "operator*(TinyVector, Matrix): Shape mismatch.");
1783 
1784  TinyVector<T, N> res;
1785  for(MultiArrayIndex i = 0; i < N; ++i)
1786  res[i] = dot(a, TinyVectorView<T, N>(&b(0,i)));
1787  return res;
1788 }
1789 
1790  /** perform matrix multiplication of matrices \a a and \a b.
1791  \a a and \a b must have matching shapes.
1792  The result is returned as a temporary matrix.
1793 
1794  <b>\#include</b> <vigra/matrix.hxx> or<br>
1795  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1796  Namespace: vigra::linalg
1797  */
1798 template <class T, class C1, class C2>
1799 inline TemporaryMatrix<T>
1801 {
1802  TemporaryMatrix<T> ret(rowCount(a), columnCount(b));
1803  mmul(a, b, ret);
1804  return ret;
1805 }
1806 
1807  /** divide matrix \a a by scalar \a b.
1808  The result is written into \a r. \a a and \a r must have the same shape.
1809 
1810  <b>\#include</b> <vigra/matrix.hxx> or<br>
1811  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1812  Namespace: vigra::linalg
1813  */
1814 template <class T, class C1, class C2>
1816 {
1817  const MultiArrayIndex rows = rowCount(a);
1818  const MultiArrayIndex cols = columnCount(a);
1819  vigra_precondition(rows == rowCount(r) && cols == columnCount(r),
1820  "sdiv(): Matrix sizes must agree.");
1821 
1822  for(MultiArrayIndex i = 0; i < cols; ++i)
1823  for(MultiArrayIndex j = 0; j < rows; ++j)
1824  r(j, i) = a(j, i) / b;
1825 }
1826 
1827  /** divide two matrices \a a and \a b pointwise.
1828  The result is written into \a r. All three matrices must have the same shape.
1829 
1830  <b>\#include</b> <vigra/matrix.hxx> or<br>
1831  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1832  Namespace: vigra::linalg
1833  */
1834 template <class T, class C1, class C2, class C3>
1837 {
1838  const MultiArrayIndex rrows = rowCount(r);
1839  const MultiArrayIndex rcols = columnCount(r);
1840  vigra_precondition(rrows == rowCount(a) && rcols == columnCount(a) &&
1841  rrows == rowCount(b) && rcols == columnCount(b),
1842  "pdiv(): Matrix shapes must agree.");
1843 
1844  for(MultiArrayIndex i = 0; i < rcols; ++i) {
1845  for(MultiArrayIndex j = 0; j < rrows; ++j) {
1846  r(j, i) = a(j, i) / b(j, i);
1847  }
1848  }
1849 }
1850 
1851  /** divide matrices \a a and \a b pointwise.
1852  \a a and \a b must have matching shapes.
1853  The result is returned as a temporary matrix.
1854 
1855  <b>\#include</b> <vigra/matrix.hxx> or<br>
1856  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1857  Namespace: vigra::linalg
1858  */
1859 template <class T, class C1, class C2>
1860 inline TemporaryMatrix<T>
1862 {
1863  TemporaryMatrix<T> ret(a.shape());
1864  pdiv(a, b, ret);
1865  return ret;
1866 }
1867 
1868  /** divide matrices \a a and \a b pointwise.
1869  \a a and \a b must have matching shapes.
1870  The result is returned as a temporary matrix.
1871 
1872  Usage:
1873 
1874  \code
1875  Matrix<double> a(m,n), b(m,n);
1876 
1877  Matrix<double> c = a / pointWise(b);
1878  // is equivalent to
1879  // Matrix<double> c = pdiv(a, b);
1880  \endcode
1881 
1882  <b>\#include</b> <vigra/matrix.hxx> or<br>
1883  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1884  Namespace: vigra::linalg
1885  */
1886 template <class T, class C, class U>
1887 inline TemporaryMatrix<T>
1888 operator/(const MultiArrayView<2, T, C> &a, PointWise<U> b)
1889 {
1890  return pdiv(a, b.t);
1891 }
1892 
1893  /** divide matrix \a a by scalar \a b.
1894  The result is returned as a temporary matrix.
1895 
1896  <b>\#include</b> <vigra/matrix.hxx> or<br>
1897  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1898  Namespace: vigra::linalg
1899  */
1900 template <class T, class C>
1901 inline TemporaryMatrix<T>
1903 {
1904  return TemporaryMatrix<T>(a) /= b;
1905 }
1906 
1907 template <class T>
1908 inline TemporaryMatrix<T>
1909 operator/(const TemporaryMatrix<T> &a, T b)
1910 {
1911  return const_cast<TemporaryMatrix<T> &>(a) /= b;
1912 }
1913 
1914  /** Create a matrix whose elements are the quotients between scalar \a a and
1915  matrix \a b. The result is returned as a temporary matrix.
1916 
1917  <b>\#include</b> <vigra/matrix.hxx> or<br>
1918  <b>\#include</b> <vigra/linear_algebra.hxx><br>
1919  Namespace: vigra::linalg
1920  */
1921 template <class T, class C>
1922 inline TemporaryMatrix<T>
1924 {
1925  return TemporaryMatrix<T>(b.shape(), a) / pointWise(b);
1926 }
1927 
1928 using vigra::argMin;
1929 using vigra::argMinIf;
1930 using vigra::argMax;
1931 using vigra::argMaxIf;
1932 
1933  /** \brief Find the index of the minimum element in a matrix.
1934 
1935  The function returns the index in column-major scan-order sense,
1936  i.e. according to the order used by <tt>MultiArrayView::operator[]</tt>.
1937  If the matrix is actually a vector, this is just the row or columns index.
1938  In case of a truly 2-dimensional matrix, the index can be converted to an
1939  index pair by calling <tt>MultiArrayView::scanOrderIndexToCoordinate()</tt>
1940 
1941  <b>Required Interface:</b>
1942 
1943  \code
1944  bool f = a[0] < NumericTraits<T>::max();
1945  \endcode
1946 
1947  <b>\#include</b> <vigra/matrix.hxx><br>
1948  Namespace: vigra
1949  */
1950 template <class T, class C>
1952 {
1953  T vopt = NumericTraits<T>::max();
1954  int best = -1;
1955  for(int k=0; k < a.size(); ++k)
1956  {
1957  if(a[k] < vopt)
1958  {
1959  vopt = a[k];
1960  best = k;
1961  }
1962  }
1963  return best;
1964 }
1965 
1966  /** \brief Find the index of the maximum element in a matrix.
1967 
1968  The function returns the index in column-major scan-order sense,
1969  i.e. according to the order used by <tt>MultiArrayView::operator[]</tt>.
1970  If the matrix is actually a vector, this is just the row or columns index.
1971  In case of a truly 2-dimensional matrix, the index can be converted to an
1972  index pair by calling <tt>MultiArrayView::scanOrderIndexToCoordinate()</tt>
1973 
1974  <b>Required Interface:</b>
1975 
1976  \code
1977  bool f = NumericTraits<T>::min() < a[0];
1978  \endcode
1979 
1980  <b>\#include</b> <vigra/matrix.hxx><br>
1981  Namespace: vigra
1982  */
1983 template <class T, class C>
1985 {
1986  T vopt = NumericTraits<T>::min();
1987  int best = -1;
1988  for(int k=0; k < a.size(); ++k)
1989  {
1990  if(vopt < a[k])
1991  {
1992  vopt = a[k];
1993  best = k;
1994  }
1995  }
1996  return best;
1997 }
1998 
1999  /** \brief Find the index of the minimum element in a matrix subject to a condition.
2000 
2001  The function returns <tt>-1</tt> if no element conforms to \a condition.
2002  Otherwise, the index of the maximum element is returned in column-major scan-order sense,
2003  i.e. according to the order used by <tt>MultiArrayView::operator[]</tt>.
2004  If the matrix is actually a vector, this is just the row or columns index.
2005  In case of a truly 2-dimensional matrix, the index can be converted to an
2006  index pair by calling <tt>MultiArrayView::scanOrderIndexToCoordinate()</tt>
2007 
2008  <b>Required Interface:</b>
2009 
2010  \code
2011  bool c = condition(a[0]);
2012  bool f = a[0] < NumericTraits<T>::max();
2013  \endcode
2014 
2015  <b>\#include</b> <vigra/matrix.hxx><br>
2016  Namespace: vigra
2017  */
2018 template <class T, class C, class UnaryFunctor>
2019 int argMinIf(MultiArrayView<2, T, C> const & a, UnaryFunctor condition)
2020 {
2021  T vopt = NumericTraits<T>::max();
2022  int best = -1;
2023  for(int k=0; k < a.size(); ++k)
2024  {
2025  if(condition(a[k]) && a[k] < vopt)
2026  {
2027  vopt = a[k];
2028  best = k;
2029  }
2030  }
2031  return best;
2032 }
2033 
2034  /** \brief Find the index of the maximum element in a matrix subject to a condition.
2035 
2036  The function returns <tt>-1</tt> if no element conforms to \a condition.
2037  Otherwise, the index of the maximum element is returned in column-major scan-order sense,
2038  i.e. according to the order used by <tt>MultiArrayView::operator[]</tt>.
2039  If the matrix is actually a vector, this is just the row or columns index.
2040  In case of a truly 2-dimensional matrix, the index can be converted to an
2041  index pair by calling <tt>MultiArrayView::scanOrderIndexToCoordinate()</tt>
2042 
2043  <b>Required Interface:</b>
2044 
2045  \code
2046  bool c = condition(a[0]);
2047  bool f = NumericTraits<T>::min() < a[0];
2048  \endcode
2049 
2050  <b>\#include</b> <vigra/matrix.hxx><br>
2051  Namespace: vigra
2052  */
2053 template <class T, class C, class UnaryFunctor>
2054 int argMaxIf(MultiArrayView<2, T, C> const & a, UnaryFunctor condition)
2055 {
2056  T vopt = NumericTraits<T>::min();
2057  int best = -1;
2058  for(int k=0; k < a.size(); ++k)
2059  {
2060  if(condition(a[k]) && vopt < a[k])
2061  {
2062  vopt = a[k];
2063  best = k;
2064  }
2065  }
2066  return best;
2067 }
2068 
2069 /** Matrix point-wise power.
2070 */
2071 template <class T, class C>
2072 linalg::TemporaryMatrix<T> pow(MultiArrayView<2, T, C> const & v, T exponent)
2073 {
2074  linalg::TemporaryMatrix<T> t(v.shape());
2075  MultiArrayIndex m = rowCount(v), n = columnCount(v);
2076 
2077  for(MultiArrayIndex i = 0; i < n; ++i)
2078  for(MultiArrayIndex j = 0; j < m; ++j)
2079  t(j, i) = vigra::pow(v(j, i), exponent);
2080  return t;
2081 }
2082 
2083 template <class T>
2084 linalg::TemporaryMatrix<T> pow(linalg::TemporaryMatrix<T> const & v, T exponent)
2085 {
2086  linalg::TemporaryMatrix<T> & t = const_cast<linalg::TemporaryMatrix<T> &>(v);
2087  MultiArrayIndex m = rowCount(t), n = columnCount(t);
2088 
2089  for(MultiArrayIndex i = 0; i < n; ++i)
2090  for(MultiArrayIndex j = 0; j < m; ++j)
2091  t(j, i) = vigra::pow(t(j, i), exponent);
2092  return t;
2093 }
2094 
2095 template <class T, class C>
2096 linalg::TemporaryMatrix<T> pow(MultiArrayView<2, T, C> const & v, int exponent)
2097 {
2098  linalg::TemporaryMatrix<T> t(v.shape());
2099  MultiArrayIndex m = rowCount(v), n = columnCount(v);
2100 
2101  for(MultiArrayIndex i = 0; i < n; ++i)
2102  for(MultiArrayIndex j = 0; j < m; ++j)
2103  t(j, i) = vigra::pow(v(j, i), exponent);
2104  return t;
2105 }
2106 
2107 template <class T>
2108 linalg::TemporaryMatrix<T> pow(linalg::TemporaryMatrix<T> const & v, int exponent)
2109 {
2110  linalg::TemporaryMatrix<T> & t = const_cast<linalg::TemporaryMatrix<T> &>(v);
2111  MultiArrayIndex m = rowCount(t), n = columnCount(t);
2112 
2113  for(MultiArrayIndex i = 0; i < n; ++i)
2114  for(MultiArrayIndex j = 0; j < m; ++j)
2115  t(j, i) = vigra::pow(t(j, i), exponent);
2116  return t;
2117 }
2118 
2119 template <class C>
2120 linalg::TemporaryMatrix<int> pow(MultiArrayView<2, int, C> const & v, int exponent)
2121 {
2122  linalg::TemporaryMatrix<int> t(v.shape());
2123  MultiArrayIndex m = rowCount(v), n = columnCount(v);
2124 
2125  for(MultiArrayIndex i = 0; i < n; ++i)
2126  for(MultiArrayIndex j = 0; j < m; ++j)
2127  t(j, i) = (int)vigra::pow((double)v(j, i), exponent);
2128  return t;
2129 }
2130 
2131 inline
2132 linalg::TemporaryMatrix<int> pow(linalg::TemporaryMatrix<int> const & v, int exponent)
2133 {
2134  linalg::TemporaryMatrix<int> & t = const_cast<linalg::TemporaryMatrix<int> &>(v);
2135  MultiArrayIndex m = rowCount(t), n = columnCount(t);
2136 
2137  for(MultiArrayIndex i = 0; i < n; ++i)
2138  for(MultiArrayIndex j = 0; j < m; ++j)
2139  t(j, i) = (int)vigra::pow((double)t(j, i), exponent);
2140  return t;
2141 }
2142 
2143  /** Matrix point-wise sqrt. */
2144 template <class T, class C>
2145 linalg::TemporaryMatrix<T> sqrt(MultiArrayView<2, T, C> const & v);
2146  /** Matrix point-wise exp. */
2147 template <class T, class C>
2148 linalg::TemporaryMatrix<T> exp(MultiArrayView<2, T, C> const & v);
2149  /** Matrix point-wise log. */
2150 template <class T, class C>
2151 linalg::TemporaryMatrix<T> log(MultiArrayView<2, T, C> const & v);
2152  /** Matrix point-wise log10. */
2153 template <class T, class C>
2154 linalg::TemporaryMatrix<T> log10(MultiArrayView<2, T, C> const & v);
2155  /** Matrix point-wise sin. */
2156 template <class T, class C>
2157 linalg::TemporaryMatrix<T> sin(MultiArrayView<2, T, C> const & v);
2158  /** Matrix point-wise asin. */
2159 template <class T, class C>
2160 linalg::TemporaryMatrix<T> asin(MultiArrayView<2, T, C> const & v);
2161  /** Matrix point-wise cos. */
2162 template <class T, class C>
2163 linalg::TemporaryMatrix<T> cos(MultiArrayView<2, T, C> const & v);
2164  /** Matrix point-wise acos. */
2165 template <class T, class C>
2166 linalg::TemporaryMatrix<T> acos(MultiArrayView<2, T, C> const & v);
2167  /** Matrix point-wise tan. */
2168 template <class T, class C>
2169 linalg::TemporaryMatrix<T> tan(MultiArrayView<2, T, C> const & v);
2170  /** Matrix point-wise atan. */
2171 template <class T, class C>
2172 linalg::TemporaryMatrix<T> atan(MultiArrayView<2, T, C> const & v);
2173  /** Matrix point-wise round. */
2174 template <class T, class C>
2175 linalg::TemporaryMatrix<T> round(MultiArrayView<2, T, C> const & v);
2176  /** Matrix point-wise floor. */
2177 template <class T, class C>
2178 linalg::TemporaryMatrix<T> floor(MultiArrayView<2, T, C> const & v);
2179  /** Matrix point-wise ceil. */
2180 template <class T, class C>
2181 linalg::TemporaryMatrix<T> ceil(MultiArrayView<2, T, C> const & v);
2182  /** Matrix point-wise abs. */
2183 template <class T, class C>
2184 linalg::TemporaryMatrix<T> abs(MultiArrayView<2, T, C> const & v);
2185  /** Matrix point-wise square. */
2186 template <class T, class C>
2187 linalg::TemporaryMatrix<T> sq(MultiArrayView<2, T, C> const & v);
2188  /** Matrix point-wise sign. */
2189 template <class T, class C>
2190 linalg::TemporaryMatrix<T> sign(MultiArrayView<2, T, C> const & v);
2191 
2192 #define VIGRA_MATRIX_UNARY_FUNCTION(FUNCTION, NAMESPACE) \
2193 using NAMESPACE::FUNCTION; \
2194 template <class T, class C> \
2195 linalg::TemporaryMatrix<T> FUNCTION(MultiArrayView<2, T, C> const & v) \
2196 { \
2197  linalg::TemporaryMatrix<T> t(v.shape()); \
2198  MultiArrayIndex m = rowCount(v), n = columnCount(v); \
2199  \
2200  for(MultiArrayIndex i = 0; i < n; ++i) \
2201  for(MultiArrayIndex j = 0; j < m; ++j) \
2202  t(j, i) = NAMESPACE::FUNCTION(v(j, i)); \
2203  return t; \
2204 } \
2205  \
2206 template <class T> \
2207 linalg::TemporaryMatrix<T> FUNCTION(linalg::Matrix<T> const & v) \
2208 { \
2209  linalg::TemporaryMatrix<T> t(v.shape()); \
2210  MultiArrayIndex m = rowCount(v), n = columnCount(v); \
2211  \
2212  for(MultiArrayIndex i = 0; i < n; ++i) \
2213  for(MultiArrayIndex j = 0; j < m; ++j) \
2214  t(j, i) = NAMESPACE::FUNCTION(v(j, i)); \
2215  return t; \
2216 } \
2217  \
2218 template <class T> \
2219 linalg::TemporaryMatrix<T> FUNCTION(linalg::TemporaryMatrix<T> const & v) \
2220 { \
2221  linalg::TemporaryMatrix<T> & t = const_cast<linalg::TemporaryMatrix<T> &>(v); \
2222  MultiArrayIndex m = rowCount(t), n = columnCount(t); \
2223  \
2224  for(MultiArrayIndex i = 0; i < n; ++i) \
2225  for(MultiArrayIndex j = 0; j < m; ++j) \
2226  t(j, i) = NAMESPACE::FUNCTION(t(j, i)); \
2227  return v; \
2228 }\
2229 }\
2230 using linalg::FUNCTION;\
2231 namespace linalg {
2232 
2233 VIGRA_MATRIX_UNARY_FUNCTION(sqrt, std)
2234 VIGRA_MATRIX_UNARY_FUNCTION(exp, std)
2235 VIGRA_MATRIX_UNARY_FUNCTION(log, std)
2236 VIGRA_MATRIX_UNARY_FUNCTION(log10, std)
2237 VIGRA_MATRIX_UNARY_FUNCTION(sin, std)
2238 VIGRA_MATRIX_UNARY_FUNCTION(asin, std)
2239 VIGRA_MATRIX_UNARY_FUNCTION(cos, std)
2240 VIGRA_MATRIX_UNARY_FUNCTION(acos, std)
2241 VIGRA_MATRIX_UNARY_FUNCTION(tan, std)
2242 VIGRA_MATRIX_UNARY_FUNCTION(atan, std)
2243 VIGRA_MATRIX_UNARY_FUNCTION(round, vigra)
2244 VIGRA_MATRIX_UNARY_FUNCTION(floor, vigra)
2245 VIGRA_MATRIX_UNARY_FUNCTION(ceil, vigra)
2246 VIGRA_MATRIX_UNARY_FUNCTION(abs, vigra)
2247 VIGRA_MATRIX_UNARY_FUNCTION(sq, vigra)
2248 VIGRA_MATRIX_UNARY_FUNCTION(sign, vigra)
2249 
2250 #undef VIGRA_MATRIX_UNARY_FUNCTION
2251 
2252 //@}
2253 
2254 } // namespace linalg
2255 
2256 using linalg::RowMajor;
2257 using linalg::ColumnMajor;
2258 using linalg::Matrix;
2261 using linalg::transpose;
2262 using linalg::pointWise;
2263 using linalg::trace;
2264 using linalg::dot;
2265 using linalg::cross;
2266 using linalg::outer;
2267 using linalg::rowCount;
2268 using linalg::columnCount;
2269 using linalg::rowVector;
2270 using linalg::columnVector;
2271 using linalg::subVector;
2272 using linalg::isSymmetric;
2275 using linalg::argMin;
2276 using linalg::argMinIf;
2277 using linalg::argMax;
2278 using linalg::argMaxIf;
2279 
2280 /********************************************************/
2281 /* */
2282 /* NormTraits */
2283 /* */
2284 /********************************************************/
2285 
2286 template <class T, class ALLOC>
2287 struct NormTraits<Matrix<T, ALLOC> >
2288 : public NormTraits<MultiArray<2, T, ALLOC> >
2289 {
2290  typedef NormTraits<MultiArray<2, T, ALLOC> > BaseType;
2291  typedef Matrix<T, ALLOC> Type;
2292  typedef typename BaseType::SquaredNormType SquaredNormType;
2293  typedef typename BaseType::NormType NormType;
2294 };
2295 
2296 template <class T, class ALLOC>
2297 struct NormTraits<linalg::TemporaryMatrix<T, ALLOC> >
2298 : public NormTraits<Matrix<T, ALLOC> >
2299 {
2300  typedef NormTraits<Matrix<T, ALLOC> > BaseType;
2301  typedef linalg::TemporaryMatrix<T, ALLOC> Type;
2302  typedef typename BaseType::SquaredNormType SquaredNormType;
2303  typedef typename BaseType::NormType NormType;
2304 };
2305 
2306 } // namespace vigra
2307 
2308 namespace std {
2309 
2310 /** \addtogroup LinearAlgebraFunctions
2311  */
2312 //@{
2313 
2314  /** print a matrix \a m to the stream \a s.
2315 
2316  <b>\#include</b> <vigra/matrix.hxx> or<br>
2317  <b>\#include</b> <vigra/linear_algebra.hxx><br>
2318  Namespace: std
2319  */
2320 template <class T, class C>
2321 ostream &
2322 operator<<(ostream & s, const vigra::MultiArrayView<2, T, C> &m)
2323 {
2326  ios::fmtflags flags = s.setf(ios::right | ios::fixed, ios::adjustfield | ios::floatfield);
2327  for(vigra::MultiArrayIndex j = 0; j < rows; ++j)
2328  {
2329  for(vigra::MultiArrayIndex i = 0; i < cols; ++i)
2330  {
2331  s << m(j, i) << " ";
2332  }
2333  s << endl;
2334  }
2335  s.setf(flags);
2336  return s;
2337 }
2338 
2339 //@}
2340 
2341 } // namespace std
2342 
2343 namespace vigra {
2344 
2345 namespace linalg {
2346 
2347 namespace detail {
2348 
2349 template <class T1, class C1, class T2, class C2, class T3, class C3>
2350 void
2351 columnStatisticsImpl(MultiArrayView<2, T1, C1> const & A,
2352  MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & sumOfSquaredDifferences)
2353 {
2354  MultiArrayIndex m = rowCount(A);
2356  vigra_precondition(1 == rowCount(mean) && n == columnCount(mean) &&
2357  1 == rowCount(sumOfSquaredDifferences) && n == columnCount(sumOfSquaredDifferences),
2358  "columnStatistics(): Shape mismatch between input and output.");
2359 
2360  // West's algorithm for incremental variance computation
2361  mean.init(NumericTraits<T2>::zero());
2362  sumOfSquaredDifferences.init(NumericTraits<T3>::zero());
2363 
2364  for(MultiArrayIndex k=0; k<m; ++k)
2365  {
2366  typedef typename NumericTraits<T2>::RealPromote TmpType;
2367  Matrix<T2> t = rowVector(A, k) - mean;
2368  TmpType f = TmpType(1.0 / (k + 1.0)),
2369  f1 = TmpType(1.0 - f);
2370  mean += f*t;
2371  sumOfSquaredDifferences += f1*sq(t);
2372  }
2373 }
2374 
2375 template <class T1, class C1, class T2, class C2, class T3, class C3>
2376 void
2377 columnStatistics2PassImpl(MultiArrayView<2, T1, C1> const & A,
2378  MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & sumOfSquaredDifferences)
2379 {
2380  MultiArrayIndex m = rowCount(A);
2382  vigra_precondition(1 == rowCount(mean) && n == columnCount(mean) &&
2383  1 == rowCount(sumOfSquaredDifferences) && n == columnCount(sumOfSquaredDifferences),
2384  "columnStatistics(): Shape mismatch between input and output.");
2385 
2386  // two-pass algorithm for incremental variance computation
2387  mean.init(NumericTraits<T2>::zero());
2388  for(MultiArrayIndex k=0; k<m; ++k)
2389  {
2390  mean += rowVector(A, k);
2391  }
2392  mean /= (double)m;
2393 
2394  sumOfSquaredDifferences.init(NumericTraits<T3>::zero());
2395  for(MultiArrayIndex k=0; k<m; ++k)
2396  {
2397  sumOfSquaredDifferences += sq(rowVector(A, k) - mean);
2398  }
2399 }
2400 
2401 } // namespace detail
2402 
2403 /** \addtogroup LinearAlgebraFunctions
2404  */
2405 //@{
2406  /** Compute statistics of every column of matrix \a A.
2407 
2408  The result matrices must be row vectors with as many columns as \a A.
2409 
2410  <b> Declarations:</b>
2411 
2412  compute only the mean:
2413  \code
2414  namespace vigra { namespace linalg {
2415  template <class T1, class C1, class T2, class C2>
2416  void
2417  columnStatistics(MultiArrayView<2, T1, C1> const & A,
2418  MultiArrayView<2, T2, C2> & mean);
2419  } }
2420  \endcode
2421 
2422  compute mean and standard deviation:
2423  \code
2424  namespace vigra { namespace linalg {
2425  template <class T1, class C1, class T2, class C2, class T3, class C3>
2426  void
2427  columnStatistics(MultiArrayView<2, T1, C1> const & A,
2428  MultiArrayView<2, T2, C2> & mean,
2429  MultiArrayView<2, T3, C3> & stdDev);
2430  } }
2431  \endcode
2432 
2433  compute mean, standard deviation, and norm:
2434  \code
2435  namespace vigra { namespace linalg {
2436  template <class T1, class C1, class T2, class C2, class T3, class C3, class T4, class C4>
2437  void
2438  columnStatistics(MultiArrayView<2, T1, C1> const & A,
2439  MultiArrayView<2, T2, C2> & mean,
2440  MultiArrayView<2, T3, C3> & stdDev,
2441  MultiArrayView<2, T4, C4> & norm);
2442  } }
2443  \endcode
2444 
2445  <b> Usage:</b>
2446 
2447  <b>\#include</b> <vigra/matrix.hxx> or<br>
2448  <b>\#include</b> <vigra/linear_algebra.hxx><br>
2449  Namespaces: vigra and vigra::linalg
2450 
2451  \code
2452  Matrix A(rows, columns);
2453  .. // fill A
2454  Matrix columnMean(1, columns), columnStdDev(1, columns), columnNorm(1, columns);
2455 
2456  columnStatistics(A, columnMean, columnStdDev, columnNorm);
2457 
2458  \endcode
2459  */
2460 doxygen_overloaded_function(template <...> void columnStatistics)
2461 
2462 template <class T1, class C1, class T2, class C2>
2463 void
2464 columnStatistics(MultiArrayView<2, T1, C1> const & A,
2465  MultiArrayView<2, T2, C2> & mean)
2466 {
2467  MultiArrayIndex m = rowCount(A);
2469  vigra_precondition(1 == rowCount(mean) && n == columnCount(mean),
2470  "columnStatistics(): Shape mismatch between input and output.");
2471 
2472  mean.init(NumericTraits<T2>::zero());
2473 
2474  for(MultiArrayIndex k=0; k<m; ++k)
2475  {
2476  mean += rowVector(A, k);
2477  }
2478  mean /= T2(m);
2479 }
2480 
2481 template <class T1, class C1, class T2, class C2, class T3, class C3>
2482 void
2483 columnStatistics(MultiArrayView<2, T1, C1> const & A,
2484  MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & stdDev)
2485 {
2486  detail::columnStatisticsImpl(A, mean, stdDev);
2487 
2488  if(rowCount(A) > 1)
2489  stdDev = sqrt(stdDev / T3(rowCount(A) - 1.0));
2490 }
2491 
2492 template <class T1, class C1, class T2, class C2, class T3, class C3, class T4, class C4>
2493 void
2494 columnStatistics(MultiArrayView<2, T1, C1> const & A,
2495  MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & stdDev, MultiArrayView<2, T4, C4> & norm)
2496 {
2497  MultiArrayIndex m = rowCount(A);
2499  vigra_precondition(1 == rowCount(mean) && n == columnCount(mean) &&
2500  1 == rowCount(stdDev) && n == columnCount(stdDev) &&
2501  1 == rowCount(norm) && n == columnCount(norm),
2502  "columnStatistics(): Shape mismatch between input and output.");
2503 
2504  detail::columnStatisticsImpl(A, mean, stdDev);
2505  norm = sqrt(stdDev + T2(m) * sq(mean));
2506  stdDev = sqrt(stdDev / T3(m - 1.0));
2507 }
2508 
2509  /** Compute statistics of every row of matrix \a A.
2510 
2511  The result matrices must be column vectors with as many rows as \a A.
2512 
2513  <b> Declarations:</b>
2514 
2515  compute only the mean:
2516  \code
2517  namespace vigra { namespace linalg {
2518  template <class T1, class C1, class T2, class C2>
2519  void
2520  rowStatistics(MultiArrayView<2, T1, C1> const & A,
2521  MultiArrayView<2, T2, C2> & mean);
2522  } }
2523  \endcode
2524 
2525  compute mean and standard deviation:
2526  \code
2527  namespace vigra { namespace linalg {
2528  template <class T1, class C1, class T2, class C2, class T3, class C3>
2529  void
2530  rowStatistics(MultiArrayView<2, T1, C1> const & A,
2531  MultiArrayView<2, T2, C2> & mean,
2532  MultiArrayView<2, T3, C3> & stdDev);
2533  } }
2534  \endcode
2535 
2536  compute mean, standard deviation, and norm:
2537  \code
2538  namespace vigra { namespace linalg {
2539  template <class T1, class C1, class T2, class C2, class T3, class C3, class T4, class C4>
2540  void
2541  rowStatistics(MultiArrayView<2, T1, C1> const & A,
2542  MultiArrayView<2, T2, C2> & mean,
2543  MultiArrayView<2, T3, C3> & stdDev,
2544  MultiArrayView<2, T4, C4> & norm);
2545  } }
2546  \endcode
2547 
2548  <b> Usage:</b>
2549 
2550  <b>\#include</b> <vigra/matrix.hxx> or<br>
2551  <b>\#include</b> <vigra/linear_algebra.hxx><br>
2552  Namespaces: vigra and vigra::linalg
2553 
2554  \code
2555  Matrix A(rows, columns);
2556  .. // fill A
2557  Matrix rowMean(rows, 1), rowStdDev(rows, 1), rowNorm(rows, 1);
2558 
2559  rowStatistics(a, rowMean, rowStdDev, rowNorm);
2560 
2561  \endcode
2562  */
2563 doxygen_overloaded_function(template <...> void rowStatistics)
2564 
2565 template <class T1, class C1, class T2, class C2>
2566 void
2567 rowStatistics(MultiArrayView<2, T1, C1> const & A,
2568  MultiArrayView<2, T2, C2> & mean)
2569 {
2570  vigra_precondition(1 == columnCount(mean) && rowCount(A) == rowCount(mean),
2571  "rowStatistics(): Shape mismatch between input and output.");
2572  MultiArrayView<2, T2, StridedArrayTag> tm = transpose(mean);
2573  columnStatistics(transpose(A), tm);
2574 }
2575 
2576 template <class T1, class C1, class T2, class C2, class T3, class C3>
2577 void
2578 rowStatistics(MultiArrayView<2, T1, C1> const & A,
2579  MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & stdDev)
2580 {
2581  vigra_precondition(1 == columnCount(mean) && rowCount(A) == rowCount(mean) &&
2582  1 == columnCount(stdDev) && rowCount(A) == rowCount(stdDev),
2583  "rowStatistics(): Shape mismatch between input and output.");
2584  MultiArrayView<2, T2, StridedArrayTag> tm = transpose(mean);
2585  MultiArrayView<2, T3, StridedArrayTag> ts = transpose(stdDev);
2586  columnStatistics(transpose(A), tm, ts);
2587 }
2588 
2589 template <class T1, class C1, class T2, class C2, class T3, class C3, class T4, class C4>
2590 void
2591 rowStatistics(MultiArrayView<2, T1, C1> const & A,
2592  MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & stdDev, MultiArrayView<2, T4, C4> & norm)
2593 {
2594  vigra_precondition(1 == columnCount(mean) && rowCount(A) == rowCount(mean) &&
2595  1 == columnCount(stdDev) && rowCount(A) == rowCount(stdDev) &&
2596  1 == columnCount(norm) && rowCount(A) == rowCount(norm),
2597  "rowStatistics(): Shape mismatch between input and output.");
2598  MultiArrayView<2, T2, StridedArrayTag> tm = transpose(mean);
2599  MultiArrayView<2, T3, StridedArrayTag> ts = transpose(stdDev);
2600  MultiArrayView<2, T4, StridedArrayTag> tn = transpose(norm);
2601  columnStatistics(transpose(A), tm, ts, tn);
2602 }
2603 
2604 namespace detail {
2605 
2606 template <class T1, class C1, class U, class T2, class C2, class T3, class C3>
2607 void updateCovarianceMatrix(MultiArrayView<2, T1, C1> const & features,
2608  U & count, MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & covariance)
2609 {
2610  MultiArrayIndex n = std::max(rowCount(features), columnCount(features));
2611  vigra_precondition(std::min(rowCount(features), columnCount(features)) == 1,
2612  "updateCovarianceMatrix(): Features must be a row or column vector.");
2613  vigra_precondition(mean.shape() == features.shape(),
2614  "updateCovarianceMatrix(): Shape mismatch between feature vector and mean vector.");
2615  vigra_precondition(n == rowCount(covariance) && n == columnCount(covariance),
2616  "updateCovarianceMatrix(): Shape mismatch between feature vector and covariance matrix.");
2617 
2618  // West's algorithm for incremental covariance matrix computation
2619  Matrix<T2> t = features - mean;
2620  ++count;
2621  double f = 1.0 / count,
2622  f1 = 1.0 - f;
2623  mean += f*t;
2624 
2625  if(rowCount(features) == 1) // update column covariance from current row
2626  {
2627  for(MultiArrayIndex k=0; k<n; ++k)
2628  {
2629  covariance(k, k) += f1*sq(t(0, k));
2630  for(MultiArrayIndex l=k+1; l<n; ++l)
2631  {
2632  covariance(k, l) += f1*t(0, k)*t(0, l);
2633  covariance(l, k) = covariance(k, l);
2634  }
2635  }
2636  }
2637  else // update row covariance from current column
2638  {
2639  for(MultiArrayIndex k=0; k<n; ++k)
2640  {
2641  covariance(k, k) += f1*sq(t(k, 0));
2642  for(MultiArrayIndex l=k+1; l<n; ++l)
2643  {
2644  covariance(k, l) += f1*t(k, 0)*t(l, 0);
2645  covariance(l, k) = covariance(k, l);
2646  }
2647  }
2648  }
2649 }
2650 
2651 } // namespace detail
2652 
2653  /** \brief Compute the covariance matrix between the columns of a matrix \a features.
2654 
2655  The result matrix \a covariance must by a square matrix with as many rows and
2656  columns as the number of columns in matrix \a features.
2657 
2658  <b>\#include</b> <vigra/matrix.hxx><br>
2659  Namespace: vigra
2660  */
2661 template <class T1, class C1, class T2, class C2>
2663  MultiArrayView<2, T2, C2> & covariance)
2664 {
2665  MultiArrayIndex m = rowCount(features), n = columnCount(features);
2666  vigra_precondition(n == rowCount(covariance) && n == columnCount(covariance),
2667  "covarianceMatrixOfColumns(): Shape mismatch between feature matrix and covariance matrix.");
2668  MultiArrayIndex count = 0;
2669  Matrix<T2> means(1, n);
2670  covariance.init(NumericTraits<T2>::zero());
2671  for(MultiArrayIndex k=0; k<m; ++k)
2672  detail::updateCovarianceMatrix(rowVector(features, k), count, means, covariance);
2673  covariance /= T2(m - 1);
2674 }
2675 
2676  /** \brief Compute the covariance matrix between the columns of a matrix \a features.
2677 
2678  The result is returned as a square temporary matrix with as many rows and
2679  columns as the number of columns in matrix \a features.
2680 
2681  <b>\#include</b> <vigra/matrix.hxx><br>
2682  Namespace: vigra
2683  */
2684 template <class T, class C>
2685 TemporaryMatrix<T>
2687 {
2688  TemporaryMatrix<T> res(columnCount(features), columnCount(features));
2689  covarianceMatrixOfColumns(features, res);
2690  return res;
2691 }
2692 
2693  /** \brief Compute the covariance matrix between the rows of a matrix \a features.
2694 
2695  The result matrix \a covariance must by a square matrix with as many rows and
2696  columns as the number of rows in matrix \a features.
2697 
2698  <b>\#include</b> <vigra/matrix.hxx><br>
2699  Namespace: vigra
2700  */
2701 template <class T1, class C1, class T2, class C2>
2703  MultiArrayView<2, T2, C2> & covariance)
2704 {
2705  MultiArrayIndex m = rowCount(features), n = columnCount(features);
2706  vigra_precondition(m == rowCount(covariance) && m == columnCount(covariance),
2707  "covarianceMatrixOfRows(): Shape mismatch between feature matrix and covariance matrix.");
2708  MultiArrayIndex count = 0;
2709  Matrix<T2> means(m, 1);
2710  covariance.init(NumericTraits<T2>::zero());
2711  for(MultiArrayIndex k=0; k<n; ++k)
2712  detail::updateCovarianceMatrix(columnVector(features, k), count, means, covariance);
2713  covariance /= T2(m - 1);
2714 }
2715 
2716  /** \brief Compute the covariance matrix between the rows of a matrix \a features.
2717 
2718  The result is returned as a square temporary matrix with as many rows and
2719  columns as the number of rows in matrix \a features.
2720 
2721  <b>\#include</b> <vigra/matrix.hxx><br>
2722  Namespace: vigra
2723  */
2724 template <class T, class C>
2725 TemporaryMatrix<T>
2727 {
2728  TemporaryMatrix<T> res(rowCount(features), rowCount(features));
2729  covarianceMatrixOfRows(features, res);
2730  return res;
2731 }
2732 
2733 enum DataPreparationGoals { ZeroMean = 1, UnitVariance = 2, UnitNorm = 4, UnitSum = 8 };
2734 
2735 inline DataPreparationGoals operator|(DataPreparationGoals l, DataPreparationGoals r)
2736 {
2737  return DataPreparationGoals(int(l) | int(r));
2738 }
2739 
2740 namespace detail {
2741 
2742 template <class T, class C1, class C2, class C3, class C4>
2743 void
2744 prepareDataImpl(const MultiArrayView<2, T, C1> & A,
2745  MultiArrayView<2, T, C2> & res, MultiArrayView<2, T, C3> & offset, MultiArrayView<2, T, C4> & scaling,
2746  DataPreparationGoals goals)
2747 {
2748  MultiArrayIndex m = rowCount(A);
2750  vigra_precondition(A.shape() == res.shape() &&
2751  n == columnCount(offset) && 1 == rowCount(offset) &&
2752  n == columnCount(scaling) && 1 == rowCount(scaling),
2753  "prepareDataImpl(): Shape mismatch between input and output.");
2754 
2755  if(!goals)
2756  {
2757  res = A;
2758  offset.init(NumericTraits<T>::zero());
2759  scaling.init(NumericTraits<T>::one());
2760  return;
2761  }
2762 
2763  bool zeroMean = (goals & ZeroMean) != 0;
2764  bool unitVariance = (goals & UnitVariance) != 0;
2765  bool unitNorm = (goals & UnitNorm) != 0;
2766  bool unitSum = (goals & UnitSum) != 0;
2767 
2768  if(unitSum)
2769  {
2770  vigra_precondition(goals == UnitSum,
2771  "prepareData(): Unit sum is not compatible with any other data preparation goal.");
2772 
2773  transformMultiArray(srcMultiArrayRange(A), destMultiArrayRange(scaling), FindSum<T>());
2774 
2775  offset.init(NumericTraits<T>::zero());
2776 
2777  for(MultiArrayIndex k=0; k<n; ++k)
2778  {
2779  if(scaling(0, k) != NumericTraits<T>::zero())
2780  {
2781  scaling(0, k) = NumericTraits<T>::one() / scaling(0, k);
2782  columnVector(res, k) = columnVector(A, k) * scaling(0, k);
2783  }
2784  else
2785  {
2786  scaling(0, k) = NumericTraits<T>::one();
2787  }
2788  }
2789 
2790  return;
2791  }
2792 
2793  vigra_precondition(!(unitVariance && unitNorm),
2794  "prepareData(): Unit variance and unit norm cannot be achieved at the same time.");
2795 
2796  Matrix<T> mean(1, n), sumOfSquaredDifferences(1, n);
2797  detail::columnStatisticsImpl(A, mean, sumOfSquaredDifferences);
2798 
2799  for(MultiArrayIndex k=0; k<n; ++k)
2800  {
2801  T stdDev = std::sqrt(sumOfSquaredDifferences(0, k) / T(m-1));
2802  if(closeAtTolerance(stdDev / mean(0,k), NumericTraits<T>::zero()))
2803  stdDev = NumericTraits<T>::zero();
2804  if(zeroMean && stdDev > NumericTraits<T>::zero())
2805  {
2806  columnVector(res, k) = columnVector(A, k) - mean(0,k);
2807  offset(0, k) = mean(0, k);
2808  mean(0, k) = NumericTraits<T>::zero();
2809  }
2810  else
2811  {
2812  columnVector(res, k) = columnVector(A, k);
2813  offset(0, k) = NumericTraits<T>::zero();
2814  }
2815 
2816  T norm = mean(0,k) == NumericTraits<T>::zero()
2817  ? std::sqrt(sumOfSquaredDifferences(0, k))
2818  : std::sqrt(sumOfSquaredDifferences(0, k) + T(m) * sq(mean(0,k)));
2819  if(unitNorm && norm > NumericTraits<T>::zero())
2820  {
2821  columnVector(res, k) /= norm;
2822  scaling(0, k) = NumericTraits<T>::one() / norm;
2823  }
2824  else if(unitVariance && stdDev > NumericTraits<T>::zero())
2825  {
2826  columnVector(res, k) /= stdDev;
2827  scaling(0, k) = NumericTraits<T>::one() / stdDev;
2828  }
2829  else
2830  {
2831  scaling(0, k) = NumericTraits<T>::one();
2832  }
2833  }
2834 }
2835 
2836 } // namespace detail
2837 
2838  /** \brief Standardize the columns of a matrix according to given <tt>DataPreparationGoals</tt>.
2839 
2840  For every column of the matrix \a A, this function computes mean,
2841  standard deviation, and norm. It then applies a linear transformation to the values of
2842  the column according to these statistics and the given <tt>DataPreparationGoals</tt>.
2843  The result is returned in matrix \a res which must have the same size as \a A.
2844  Optionally, the transformation applied can also be returned in the matrices \a offset
2845  and \a scaling (see below for an example how these matrices can be used to standardize
2846  more data according to the same transformation).
2847 
2848  The following <tt>DataPreparationGoals</tt> are supported:
2849 
2850  <DL>
2851  <DT><tt>ZeroMean</tt><DD> Subtract the column mean form every column if the values in the column are not constant.
2852  Do nothing in a constant column.
2853  <DT><tt>UnitSum</tt><DD> Scale the columns so that the their sum is one if the sum was initially non-zero.
2854  Do nothing in a zero-sum column.
2855  <DT><tt>UnitVariance</tt><DD> Divide by the column standard deviation if the values in the column are not constant.
2856  Do nothing in a constant column.
2857  <DT><tt>UnitNorm</tt><DD> Divide by the column norm if it is non-zero.
2858  <DT><tt>ZeroMean | UnitVariance</tt><DD> First subtract the mean and then divide by the standard deviation, unless the
2859  column is constant (in which case the column remains unchanged).
2860  <DT><tt>ZeroMean | UnitNorm</tt><DD> If the column is non-constant, subtract the mean. Then divide by the norm
2861  of the result if the norm is non-zero.
2862  </DL>
2863 
2864  <b> Declarations:</b>
2865 
2866  Standardize the matrix and return the parameters of the linear transformation.
2867  The matrices \a offset and \a scaling must be row vectors with as many columns as \a A.
2868  \code
2869  namespace vigra { namespace linalg {
2870  template <class T, class C1, class C2, class C3, class C4>
2871  void
2872  prepareColumns(MultiArrayView<2, T, C1> const & A,
2873  MultiArrayView<2, T, C2> & res,
2874  MultiArrayView<2, T, C3> & offset,
2875  MultiArrayView<2, T, C4> & scaling,
2876  DataPreparationGoals goals = ZeroMean | UnitVariance);
2877  } }
2878  \endcode
2879 
2880  Only standardize the matrix.
2881  \code
2882  namespace vigra { namespace linalg {
2883  template <class T, class C1, class C2>
2884  void
2885  prepareColumns(MultiArrayView<2, T, C1> const & A,
2886  MultiArrayView<2, T, C2> & res,
2887  DataPreparationGoals goals = ZeroMean | UnitVariance);
2888  } }
2889  \endcode
2890 
2891  <b> Usage:</b>
2892 
2893  <b>\#include</b> <vigra/matrix.hxx> or<br>
2894  <b>\#include</b> <vigra/linear_algebra.hxx><br>
2895  Namespaces: vigra and vigra::linalg
2896 
2897  \code
2898  Matrix A(rows, columns);
2899  .. // fill A
2900  Matrix standardizedA(rows, columns), offset(1, columns), scaling(1, columns);
2901 
2902  prepareColumns(A, standardizedA, offset, scaling, ZeroMean | UnitNorm);
2903 
2904  // use offset and scaling to prepare additional data according to the same transformation
2905  Matrix newData(nrows, columns);
2906 
2907  Matrix standardizedNewData = (newData - repeatMatrix(offset, nrows, 1)) * pointWise(repeatMatrix(scaling, nrows, 1));
2908 
2909  \endcode
2910  */
2911 doxygen_overloaded_function(template <...> void prepareColumns)
2912 
2913 template <class T, class C1, class C2, class C3, class C4>
2914 inline void
2915 prepareColumns(MultiArrayView<2, T, C1> const & A,
2916  MultiArrayView<2, T, C2> & res, MultiArrayView<2, T, C3> & offset, MultiArrayView<2, T, C4> & scaling,
2917  DataPreparationGoals goals = ZeroMean | UnitVariance)
2918 {
2919  detail::prepareDataImpl(A, res, offset, scaling, goals);
2920 }
2921 
2922 template <class T, class C1, class C2>
2923 inline void
2924 prepareColumns(MultiArrayView<2, T, C1> const & A, MultiArrayView<2, T, C2> & res,
2925  DataPreparationGoals goals = ZeroMean | UnitVariance)
2926 {
2927  Matrix<T> offset(1, columnCount(A)), scaling(1, columnCount(A));
2928  detail::prepareDataImpl(A, res, offset, scaling, goals);
2929 }
2930 
2931  /** \brief Standardize the rows of a matrix according to given <tt>DataPreparationGoals</tt>.
2932 
2933  This algorithm works in the same way as \ref prepareColumns() (see there for detailed
2934  documentation), but is applied to the rows of the matrix \a A instead. Accordingly, the
2935  matrices holding the parameters of the linear transformation must be column vectors
2936  with as many rows as \a A.
2937 
2938  <b> Declarations:</b>
2939 
2940  Standardize the matrix and return the parameters of the linear transformation.
2941  The matrices \a offset and \a scaling must be column vectors
2942  with as many rows as \a A.
2943 
2944  \code
2945  namespace vigra { namespace linalg {
2946  template <class T, class C1, class C2, class C3, class C4>
2947  void
2948  prepareRows(MultiArrayView<2, T, C1> const & A,
2949  MultiArrayView<2, T, C2> & res,
2950  MultiArrayView<2, T, C3> & offset,
2951  MultiArrayView<2, T, C4> & scaling,
2952  DataPreparationGoals goals = ZeroMean | UnitVariance)´;
2953  } }
2954  \endcode
2955 
2956  Only standardize the matrix.
2957  \code
2958  namespace vigra { namespace linalg {
2959  template <class T, class C1, class C2>
2960  void
2961  prepareRows(MultiArrayView<2, T, C1> const & A,
2962  MultiArrayView<2, T, C2> & res,
2963  DataPreparationGoals goals = ZeroMean | UnitVariance);
2964  } }
2965  \endcode
2966 
2967  <b> Usage:</b>
2968 
2969  <b>\#include</b> <vigra/matrix.hxx> or<br>
2970  <b>\#include</b> <vigra/linear_algebra.hxx><br>
2971  Namespaces: vigra and vigra::linalg
2972 
2973  \code
2974  Matrix A(rows, columns);
2975  .. // fill A
2976  Matrix standardizedA(rows, columns), offset(rows, 1), scaling(rows, 1);
2977 
2978  prepareRows(A, standardizedA, offset, scaling, ZeroMean | UnitNorm);
2979 
2980  // use offset and scaling to prepare additional data according to the same transformation
2981  Matrix newData(rows, ncolumns);
2982 
2983  Matrix standardizedNewData = (newData - repeatMatrix(offset, 1, ncolumns)) * pointWise(repeatMatrix(scaling, 1, ncolumns));
2984 
2985  \endcode
2986 */
2987 doxygen_overloaded_function(template <...> void prepareRows)
2988 
2989 template <class T, class C1, class C2, class C3, class C4>
2990 inline void
2991 prepareRows(MultiArrayView<2, T, C1> const & A,
2992  MultiArrayView<2, T, C2> & res, MultiArrayView<2, T, C3> & offset, MultiArrayView<2, T, C4> & scaling,
2993  DataPreparationGoals goals = ZeroMean | UnitVariance)
2994 {
2995  MultiArrayView<2, T, StridedArrayTag> tr = transpose(res), to = transpose(offset), ts = transpose(scaling);
2996  detail::prepareDataImpl(transpose(A), tr, to, ts, goals);
2997 }
2998 
2999 template <class T, class C1, class C2>
3000 inline void
3001 prepareRows(MultiArrayView<2, T, C1> const & A, MultiArrayView<2, T, C2> & res,
3002  DataPreparationGoals goals = ZeroMean | UnitVariance)
3003 {
3004  MultiArrayView<2, T, StridedArrayTag> tr = transpose(res);
3005  Matrix<T> offset(1, rowCount(A)), scaling(1, rowCount(A));
3006  detail::prepareDataImpl(transpose(A), tr, offset, scaling, goals);
3007 }
3008 
3009 //@}
3010 
3011 } // namespace linalg
3012 
3015 using linalg::rowStatistics;
3016 using linalg::prepareRows;
3017 using linalg::ZeroMean;
3018 using linalg::UnitVariance;
3019 using linalg::UnitNorm;
3020 using linalg::UnitSum;
3021 
3022 } // namespace vigra
3023 
3024 
3025 
3026 #endif // VIGRA_MATRIX_HXX
Matrix(difference_type_1 rows, difference_type_1 columns, const_reference init, allocator_type const &alloc=allocator_type())
Definition: matrix.hxx:185
Matrix & operator-=(MultiArrayView< 2, U, C > const &other)
Definition: matrix.hxx:484
Matrix & operator+=(T other)
Definition: matrix.hxx:510
Matrix & operator-=(T other)
Definition: matrix.hxx:518
linalg::TemporaryMatrix< T > acos(MultiArrayView< 2, T, C > const &v)
difference_type_1 rowCount() const
Definition: matrix.hxx:365
NormTraits< Matrix >::SquaredNormType squaredNorm() const
void covarianceMatrixOfColumns(MultiArrayView< 2, T1, C1 > const &features, MultiArrayView< 2, T2, C2 > &covariance)
Compute the covariance matrix between the columns of a matrix features.
Definition: matrix.hxx:2662
MultiArrayView< 2, T, C > columnVector(MultiArrayView< 2, T, C > const &m, MultiArrayIndex d)
Definition: matrix.hxx:725
linalg::TemporaryMatrix< T > ceil(MultiArrayView< 2, T, C > const &v)
view_type::pointer pointer
Definition: multi_array.hxx:2385
MultiArrayShape< actual_dimension >::type difference_type
Definition: multi_array.hxx:690
Matrix(difference_type_1 rows, difference_type_1 columns, const_pointer init, RawArrayMemoryLayout layout=RowMajor, allocator_type const &alloc=allocator_type())
Definition: matrix.hxx:219
MultiArrayIndex rowCount(const MultiArrayView< 2, T, C > &x)
Definition: matrix.hxx:669
void transpose(const MultiArrayView< 2, T, C1 > &v, MultiArrayView< 2, T, C2 > &r)
Definition: matrix.hxx:963
MultiArray & operator/=(const MultiArrayView< N, U, StrideTag > &rhs)
Definition: multi_array.hxx:2627
linalg::TemporaryMatrix< T > sqrt(MultiArrayView< 2, T, C > const &v)
Matrix & operator=(const Matrix &rhs)
Definition: matrix.hxx:270
void cross(const MultiArrayView< 1, T, C1 > &x, const MultiArrayView< 1, T, C2 > &y, MultiArrayView< 1, T, C3 > &r)
Definition: matrix.hxx:1402
linalg::TemporaryMatrix< T > sin(MultiArrayView< 2, T, C > const &v)
difference_type_1 columnCount() const
Definition: matrix.hxx:372
int argMin(MultiArrayView< 2, T, C > const &a)
Find the index of the minimum element in a matrix.
Definition: matrix.hxx:1951
Definition: matrix.hxx:121
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
difference_type_1 elementCount() const
Definition: matrix.hxx:379
linalg::TemporaryMatrix< T > sign(MultiArrayView< 2, T, C > const &v)
Matrix & init(const U &init)
Definition: matrix.hxx:315
TemporaryMatrix< T > operator+(const MultiArrayView< 2, T, C1 > &a, const MultiArrayView< 2, T, C2 > &b)
Definition: matrix.hxx:1134
Definition: array_vector.hxx:903
Matrix & operator/=(T other)
Definition: matrix.hxx:534
difference_type_1 elementCount() const
Definition: multi_array.hxx:1533
linalg::TemporaryMatrix< T > round(MultiArrayView< 2, T, C > const &v)
int argMinIf(MultiArrayView< 2, T, C > const &a, UnaryFunctor condition)
Find the index of the minimum element in a matrix subject to a condition.
Definition: matrix.hxx:2019
TemporaryMatrix< T > joinVertically(const MultiArrayView< 2, T, C1 > &a, const MultiArrayView< 2, T, C2 > &b)
Definition: matrix.hxx:1008
Main MultiArray class containing the memory management.
Definition: multi_array.hxx:2357
void repeatMatrix(MultiArrayView< 2, T, C1 > const &v, MultiArrayView< 2, T, C2 > &r, unsigned int verticalCount, unsigned int horizontalCount)
Definition: matrix.hxx:1061
TemporaryMatrix< T > ones(MultiArrayIndex rows, MultiArrayIndex cols)
Definition: matrix.hxx:887
void reshape(const difference_type &shape)
Definition: multi_array.hxx:2738
linalg::TemporaryMatrix< T > asin(MultiArrayView< 2, T, C > const &v)
bool isSymmetric(const MultiArrayView< 2, T, C > &v)
Definition: matrix.hxx:777
Matrix & operator*=(MultiArrayView< 2, U, C > const &other)
Definition: matrix.hxx:493
Find the average pixel value in an image or ROI.
Definition: inspectimage.hxx:1248
TemporaryMatrix< T > operator/(const MultiArrayView< 2, T, C > &a, PointWise< U > b)
Definition: matrix.hxx:1888
bool isSymmetric() const
Definition: matrix.hxx:386
void covarianceMatrixOfRows(MultiArrayView< 2, T1, C1 > const &features, MultiArrayView< 2, T2, C2 > &covariance)
Compute the covariance matrix between the rows of a matrix features.
Definition: matrix.hxx:2702
int argMaxIf(MultiArrayView< 2, T, C > const &a, UnaryFunctor condition)
Find the index of the maximum element in a matrix subject to a condition.
Definition: matrix.hxx:2054
void add(const MultiArrayView< 2, T, C1 > &a, const MultiArrayView< 2, T, C2 > &b, MultiArrayView< 2, T, C3 > &r)
Definition: matrix.hxx:1108
MultiArray & operator+=(const MultiArrayView< N, U, StrideTag > &rhs)
Definition: multi_array.hxx:2580
std::ptrdiff_t MultiArrayIndex
Definition: multi_shape.hxx:55
Find the sum of the pixel values in an image or ROI.
Definition: inspectimage.hxx:1143
void mmul(const MultiArrayView< 2, T, C1 > &a, const MultiArrayView< 2, T, C2 > &b, MultiArrayView< 2, T, C3 > &r)
Definition: matrix.hxx:1603
Definition: accessor.hxx:43
Matrix(difference_type_1 rows, difference_type_1 columns, ALLOC const &alloc=allocator_type())
Definition: matrix.hxx:165
Matrix & operator/=(MultiArrayView< 2, U, C > const &other)
Definition: matrix.hxx:502
MultiArrayView< N, T, StridedArrayTag > transpose() const
Definition: multi_array.hxx:1470
view_type::difference_type difference_type
Definition: multi_array.hxx:2405
FFTWComplex< R > & operator-=(FFTWComplex< R > &a, const FFTWComplex< R > &b)
subtract-assignment
Definition: fftw3.hxx:867
MultiArray & operator-=(const MultiArrayView< N, U, StrideTag > &rhs)
Definition: multi_array.hxx:2596
Matrix(ALLOC const &alloc)
Definition: matrix.hxx:146
NumericTraits< T >::Promote trace(MultiArrayView< 2, T, C > const &m)
Definition: matrix.hxx:799
void prepareRows(...)
Standardize the rows of a matrix according to given DataPreparationGoals.
MultiArray & operator=(const MultiArray &rhs)
Definition: multi_array.hxx:2546
MultiArrayView< N-M, T, StridedArrayTag > bindInner(const TinyVector< Index, M > &d) const
Definition: multi_array.hxx:2090
void sdiv(const MultiArrayView< 2, T, C1 > &a, T b, MultiArrayView< 2, T, C2 > &r)
Definition: matrix.hxx:1815
view_type rowVector(difference_type_1 d) const
Definition: matrix.hxx:351
difference_type_1 size() const
Definition: multi_array.hxx:1544
TemporaryMatrix< T > mean() const
Definition: matrix.hxx:416
void columnStatistics(...)
FFTWComplex< R > & operator+=(FFTWComplex< R > &a, const FFTWComplex< R > &b)
add-assignment
Definition: fftw3.hxx:859
view_type::reference reference
Definition: multi_array.hxx:2393
TemporaryMatrix< T > operator-(const MultiArrayView< 2, T, C1 > &a, const MultiArrayView< 2, T, C2 > &b)
Definition: matrix.hxx:1236
linalg::TemporaryMatrix< T > log10(MultiArrayView< 2, T, C > const &v)
void sub(const MultiArrayView< 2, T, C1 > &a, const MultiArrayView< 2, T, C2 > &b, MultiArrayView< 2, T, C3 > &r)
Definition: matrix.hxx:1210
linalg::TemporaryMatrix< T > pow(MultiArrayView< 2, T, C > const &v, T exponent)
Definition: matrix.hxx:2072
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector's elements
Definition: tinyvector.hxx:1871
TemporaryMatrix< T > joinHorizontally(const MultiArrayView< 2, T, C1 > &a, const MultiArrayView< 2, T, C2 > &b)
Definition: matrix.hxx:1034
MultiArray & init(const U &init)
Definition: multi_array.hxx:2728
TinyVector< V, SIZE > pow(TinyVectorBase< V, SIZE, D1, D2 > const &v, E exponent)
Definition: tinyvector.hxx:1834
linalg::TemporaryMatrix< T > sq(MultiArrayView< 2, T, C > const &v)
Matrix(const difference_type &shape, const_reference init, allocator_type const &alloc=allocator_type())
Definition: matrix.hxx:175
Matrix< T, ALLLOC >::NormType norm(const Matrix< T, ALLLOC > &a)
Matrix(const MultiArrayView< 2, U, C > &rhs)
Definition: matrix.hxx:261
Matrix< T, ALLLOC >::SquaredNormType squaredNorm(const Matrix< T, ALLLOC > &a)
allocator_type const & allocator() const
Definition: multi_array.hxx:2787
Matrix & operator=(const MultiArrayView< 2, U, C > &rhs)
Definition: matrix.hxx:298
MultiArrayView< 2, T, C > subVector(MultiArrayView< 2, T, C > const &m, int first, int end)
Definition: matrix.hxx:759
FFTWComplex< R > & operator*=(FFTWComplex< R > &a, const FFTWComplex< R > &b)
multiply-assignment
Definition: fftw3.hxx:875
Wrapper for fixed size vectors.
Definition: tinyvector.hxx:612
void reshape(difference_type_1 rows, difference_type_1 columns)
Definition: matrix.hxx:323
view_type::const_pointer const_pointer
Definition: multi_array.hxx:2389
view_type::value_type value_type
Definition: multi_array.hxx:2381
Matrix & operator=(const TemporaryMatrix< T, ALLOC > &rhs)
Definition: matrix.hxx:282
Matrix & operator*=(T other)
Definition: matrix.hxx:526
Class for fixed size vectors.This class contains an array of size SIZE of the specified VALUETYPE...
Definition: accessor.hxx:940
void diagonalMatrix(MultiArrayView< 2, T, C1 > const &v, MultiArrayView< 2, T, C2 > &r)
Definition: matrix.hxx:913
int argMax(MultiArrayView< 2, T, C > const &a)
Find the index of the maximum element in a matrix.
Definition: matrix.hxx:1984
NormTraits< T >::SquaredNormType dot(const MultiArrayView< 2, T, C1 > &x, const MultiArrayView< 2, T, C2 > &y)
Definition: matrix.hxx:1340
void outer(const MultiArrayView< 2, T, C1 > &x, const MultiArrayView< 2, T, C2 > &y, MultiArrayView< 2, T, C3 > &r)
Definition: matrix.hxx:1457
MultiArrayView< 2, vluae_type, StridedArrayTag > transpose() const
void reshape(difference_type_1 rows, difference_type_1 columns, const_reference init)
Definition: matrix.hxx:330
bool closeAtTolerance(T1 l, T2 r, typename PromoteTraits< T1, T2 >::Promote epsilon)
Tolerance based floating-point comparison.
Definition: mathutil.hxx:1600
linalg::TemporaryMatrix< T > log(MultiArrayView< 2, T, C > const &v)
TemporaryMatrix< T > mean(difference_type_1 d) const
Definition: matrix.hxx:427
void identityMatrix(MultiArrayView< 2, T, C > &r)
Definition: matrix.hxx:843
TemporaryMatrix< T > sum(difference_type_1 d) const
Definition: matrix.hxx:404
MultiArrayView< 2, T, C > rowVector(MultiArrayView< 2, T, C > const &m, MultiArrayIndex d)
Definition: matrix.hxx:695
void prepareColumns(...)
Standardize the columns of a matrix according to given DataPreparationGoals.
view_type::difference_type_1 difference_type_1
Definition: multi_array.hxx:2409
MultiArrayIndex columnCount(const MultiArrayView< 2, T, C > &x)
Definition: matrix.hxx:682
void pmul(const MultiArrayView< 2, T, C1 > &a, const MultiArrayView< 2, T, C2 > &b, MultiArrayView< 2, T, C3 > &r)
Definition: matrix.hxx:1648
Base class for fixed size vectors.
Definition: tinyvector.hxx:81
Matrix(const difference_type &shape, const_pointer init, RawArrayMemoryLayout layout=RowMajor, allocator_type const &alloc=allocator_type())
Definition: matrix.hxx:197
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:655
void transformMultiArray(...)
Transform a multi-dimensional array with a unary function or functor.
linalg::TemporaryMatrix< T > atan(MultiArrayView< 2, T, C > const &v)
void reshape(difference_type const &shape, const_reference init)
Definition: matrix.hxx:344
void smul(const MultiArrayView< 2, T, C1 > &a, T b, MultiArrayView< 2, T, C2 > &r)
Definition: matrix.hxx:1570
linalg::TemporaryMatrix< T > tan(MultiArrayView< 2, T, C > const &v)
MultiArrayView & init(const U &init)
Definition: multi_array.hxx:1150
void reshape(difference_type const &shape)
Definition: matrix.hxx:337
FFTWComplex< R > & operator/=(FFTWComplex< R > &a, const FFTWComplex< R > &b)
divide-assignment
Definition: fftw3.hxx:884
Matrix(const difference_type &shape, ALLOC const &alloc=allocator_type())
Definition: matrix.hxx:155
MultiArrayView subarray(difference_type p, difference_type q) const
Definition: multi_array.hxx:1431
linalg::TemporaryMatrix< T > abs(MultiArrayView< 2, T, C > const &v)
Matrix(const Matrix &rhs)
Definition: matrix.hxx:237
linalg::TemporaryMatrix< T > cos(MultiArrayView< 2, T, C > const &v)
Matrix(const TemporaryMatrix< T, ALLOC > &rhs)
Definition: matrix.hxx:251
void swap(MultiArray &other)
Matrix()
Definition: matrix.hxx:141
view_type columnVector(difference_type_1 d) const
Definition: matrix.hxx:358
NormTraits< Matrix >::NormType norm() const
TemporaryMatrix< T > sum() const
Definition: matrix.hxx:393
void rowStatistics(...)
value_type & operator()(difference_type_1 row, difference_type_1 column)
MultiArrayView< N-M, T, StrideTag > bindOuter(const TinyVector< Index, M > &d) const
Definition: multi_array.hxx:2067
void pdiv(const MultiArrayView< 2, T, C1 > &a, const MultiArrayView< 2, T, C2 > &b, MultiArrayView< 2, T, C3 > &r)
Definition: matrix.hxx:1835
Matrix & operator+=(MultiArrayView< 2, U, C > const &other)
Definition: matrix.hxx:475
linalg::TemporaryMatrix< T > floor(MultiArrayView< 2, T, C > const &v)
view_type::const_reference const_reference
Definition: multi_array.hxx:2397
Matrix & operator=(value_type const &v)
Definition: matrix.hxx:307
TemporaryMatrix< T > operator*(const MultiArrayView< 2, T, C > &a, PointWise< U > b)
Definition: matrix.hxx:1701
MultiArray & operator*=(const MultiArrayView< N, U, StrideTag > &rhs)
Definition: multi_array.hxx:2611

© 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.10.0