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

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

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