37 #ifndef VIGRA_MATRIX_HXX
38 #define VIGRA_MATRIX_HXX
43 #include "multi_array.hxx"
44 #include "mathutil.hxx"
45 #include "numerictraits.hxx"
46 #include "multi_pointoperators.hxx"
65 template <
class T,
class C>
67 rowCount(
const MultiArrayView<2, T, C> &x);
69 template <
class T,
class C>
73 template <
class T,
class C>
74 inline MultiArrayView <2, T, C>
77 template <
class T,
class C>
78 inline MultiArrayView <2, T, C>
81 template <
class T,
class ALLOC = std::allocator<T> >
82 class TemporaryMatrix;
84 template <
class T,
class C1,
class C2>
85 void transpose(
const MultiArrayView<2, T, C1> &v, MultiArrayView<2, T, C2> &r);
87 template <
class T,
class C>
90 enum RawArrayMemoryLayout { RowMajor, ColumnMajor };
120 template <
class T,
class ALLOC = std::allocator<T> >
128 typedef TemporaryMatrix<T, ALLOC> temp_type;
137 typedef ALLOC allocator_type;
156 ALLOC
const & alloc = allocator_type())
165 Matrix(difference_type_1 rows, difference_type_1 columns,
166 ALLOC
const & alloc = allocator_type())
167 :
BaseType(difference_type(rows, columns), alloc)
176 allocator_type
const & alloc = allocator_type())
185 Matrix(difference_type_1 rows, difference_type_1 columns, const_reference init,
186 allocator_type
const & alloc = allocator_type())
187 :
BaseType(difference_type(rows, columns), init, alloc)
197 Matrix(
const difference_type &shape, const_pointer init, RawArrayMemoryLayout layout = RowMajor,
198 allocator_type
const & alloc = allocator_type())
201 if(layout == RowMajor)
203 difference_type trans(shape[1], shape[0]);
219 Matrix(difference_type_1 rows, difference_type_1 columns, const_pointer init, RawArrayMemoryLayout layout = RowMajor,
220 allocator_type
const & alloc = allocator_type())
221 :
BaseType(difference_type(rows, columns), alloc)
223 if(layout == RowMajor)
225 difference_type trans(columns, rows);
251 Matrix(
const TemporaryMatrix<T, ALLOC> &rhs)
254 this->
swap(
const_cast<TemporaryMatrix<T, ALLOC> &
>(rhs));
260 template<
class U,
class C>
284 if(this->
shape() == rhs.shape())
287 this->
swap(
const_cast<TemporaryMatrix<T, ALLOC> &
>(rhs));
297 template <
class U,
class C>
315 void reshape(difference_type_1 rows, difference_type_1 columns)
322 void reshape(difference_type_1 rows, difference_type_1 columns, const_reference init)
336 void reshape(difference_type
const & shape, const_reference init)
385 TemporaryMatrix<T>
sum()
const
387 TemporaryMatrix<T> result(1, 1);
389 destMultiArrayRange(result),
396 TemporaryMatrix<T>
sum(difference_type_1 d)
const
399 TemporaryMatrix<T> result(shape);
401 destMultiArrayRange(result),
408 TemporaryMatrix<T>
mean()
const
410 TemporaryMatrix<T> result(1, 1);
412 destMultiArrayRange(result),
419 TemporaryMatrix<T>
mean(difference_type_1 d)
const
422 TemporaryMatrix<T> result(shape);
424 destMultiArrayRange(result),
437 value_type &
operator()(difference_type_1 row, difference_type_1 column);
443 value_type
operator()(difference_type_1 row, difference_type_1 column)
const;
447 typename NormTraits<Matrix>::SquaredNormType
squaredNorm()
const;
451 typename NormTraits<Matrix>::NormType
norm()
const;
466 template <
class U,
class C>
475 template <
class U,
class C>
484 template <
class U,
class C>
493 template <
class U,
class C>
537 template <
class T,
class ALLOC>
538 class TemporaryMatrix
539 :
public Matrix<T, ALLOC>
541 typedef Matrix<T, ALLOC> BaseType;
543 typedef Matrix<T, ALLOC> matrix_type;
544 typedef TemporaryMatrix<T, ALLOC> temp_type;
546 typedef typename BaseType::value_type value_type;
547 typedef typename BaseType::pointer pointer;
548 typedef typename BaseType::const_pointer const_pointer;
549 typedef typename BaseType::reference reference;
550 typedef typename BaseType::const_reference const_reference;
551 typedef typename BaseType::difference_type difference_type;
552 typedef typename BaseType::difference_type_1 difference_type_1;
553 typedef ALLOC allocator_type;
555 TemporaryMatrix(difference_type
const & shape)
556 : BaseType(shape, ALLOC())
559 TemporaryMatrix(difference_type
const & shape, const_reference init)
560 : BaseType(shape, init, ALLOC())
563 TemporaryMatrix(difference_type_1 rows, difference_type_1 columns)
564 : BaseType(rows, columns, ALLOC())
567 TemporaryMatrix(difference_type_1 rows, difference_type_1 columns, const_reference init)
568 : BaseType(rows, columns, init, ALLOC())
571 template<
class U,
class C>
572 TemporaryMatrix(
const MultiArrayView<2, U, C> &rhs)
576 TemporaryMatrix(
const TemporaryMatrix &rhs)
579 this->
swap(const_cast<TemporaryMatrix &>(rhs));
583 TemporaryMatrix & init(
const U & init)
589 template <
class U,
class C>
590 TemporaryMatrix &
operator+=(MultiArrayView<2, U, C>
const & other)
596 template <
class U,
class C>
597 TemporaryMatrix &
operator-=(MultiArrayView<2, U, C>
const & other)
603 template <
class U,
class C>
604 TemporaryMatrix &
operator*=(MultiArrayView<2, U, C>
const & other)
610 template <
class U,
class C>
611 TemporaryMatrix &
operator/=(MultiArrayView<2, U, C>
const & other)
642 TemporaryMatrix &
operator=(
const TemporaryMatrix &rhs);
659 template <
class T,
class C>
672 template <
class T,
class C>
685 template <
class T,
class C>
701 template <
class T,
class C>
706 return m.
subarray(first, Shape(first[0]+1, end));
715 template <
class T,
class C>
730 template <
class T,
class C>
735 return m.
subarray(first, Shape(end, first[1]+1));
749 template <
class T,
class C>
755 return m.
subarray(Shape(first, 0), Shape(end, 1));
756 vigra_precondition(
rowCount(m) == 1,
757 "linalg::subVector(): Input must be a vector (1xN or Nx1).");
758 return m.
subarray(Shape(0, first), Shape(1, end));
767 template <
class T,
class C>
777 if(m(j, i) != m(i, j))
789 template <
class T,
class C>
790 typename NumericTraits<T>::Promote
793 typedef typename NumericTraits<T>::Promote SumType;
796 vigra_precondition(size ==
columnCount(m),
"linalg::trace(): Matrix must be square.");
798 SumType
sum = NumericTraits<SumType>::zero();
804 #ifdef DOXYGEN // documentation only -- function is already defined in vigra/multi_array.hxx
812 template <
class T,
class ALLOC>
813 typename Matrix<T, ALLLOC>::SquaredNormType
822 template <
class T,
class ALLOC>
823 typename Matrix<T, ALLLOC>::NormType
824 norm(
const Matrix<T, ALLLOC> &a);
834 template <
class T,
class C>
839 "identityMatrix(): Matrix must be square.");
842 r(j, i) = NumericTraits<T>::zero();
843 r(i, i) = NumericTraits<T>::one();
861 TemporaryMatrix<T> ret(size, size, NumericTraits<T>::zero());
863 ret(i, i) = NumericTraits<T>::one();
881 return TemporaryMatrix<T>(rows, cols, NumericTraits<T>::one());
886 template <
class T,
class C1,
class C2>
891 "diagonalMatrix(): result must be a square matrix.");
904 template <
class T,
class C1,
class C2>
908 "diagonalMatrix(): input must be a vector.");
909 r.
init(NumericTraits<T>::zero());
932 template <
class T,
class C>
936 "diagonalMatrix(): input must be a vector.");
938 TemporaryMatrix<T> ret(size, size, NumericTraits<T>::zero());
954 template <
class T,
class C1,
class C2>
960 "transpose(): arrays must have transposed shapes.");
983 template <
class T,
class C>
998 template <
class T,
class C1,
class C2>
999 inline TemporaryMatrix<T>
1002 typedef typename TemporaryMatrix<T>::difference_type Shape;
1006 "joinVertically(): shape mismatch.");
1010 TemporaryMatrix<T> t(ma + mb, n, T());
1011 t.subarray(Shape(0,0), Shape(ma, n)) = a;
1012 t.subarray(Shape(ma,0), Shape(ma+mb, n)) = b;
1024 template <
class T,
class C1,
class C2>
1025 inline TemporaryMatrix<T>
1028 typedef typename TemporaryMatrix<T>::difference_type Shape;
1031 vigra_precondition(m ==
rowCount(b),
1032 "joinHorizontally(): shape mismatch.");
1036 TemporaryMatrix<T> t(m, na + nb, T());
1037 t.subarray(Shape(0,0), Shape(m, na)) = a;
1038 t.subarray(Shape(0, na), Shape(m, na + nb)) = b;
1052 template <
class T,
class C1,
class C2>
1054 unsigned int verticalCount,
unsigned int horizontalCount)
1060 "repeatMatrix(): Shape mismatch.");
1066 r.
subarray(Shape(k*m, l*n), Shape((k+1)*m, (l+1)*n)) = v;
1082 template <
class T,
class C>
1087 TemporaryMatrix<T> ret(verticalCount*m, horizontalCount*n);
1099 template <
class T,
class C1,
class C2,
class C3>
1107 "add(): Matrix shapes must agree.");
1111 r(j, i) = a(j, i) + b(j, i);
1124 template <
class T,
class C1,
class C2>
1125 inline TemporaryMatrix<T>
1128 return TemporaryMatrix<T>(a) += b;
1131 template <
class T,
class C>
1132 inline TemporaryMatrix<T>
1135 return const_cast<TemporaryMatrix<T> &
>(a) += b;
1138 template <
class T,
class C>
1139 inline TemporaryMatrix<T>
1140 operator+(
const MultiArrayView<2, T, C> &a,
const TemporaryMatrix<T> &b)
1142 return const_cast<TemporaryMatrix<T> &
>(b) += a;
1146 inline TemporaryMatrix<T>
1147 operator+(
const TemporaryMatrix<T> &a,
const TemporaryMatrix<T> &b)
1149 return const_cast<TemporaryMatrix<T> &
>(a) += b;
1159 template <
class T,
class C>
1160 inline TemporaryMatrix<T>
1163 return TemporaryMatrix<T>(a) += b;
1167 inline TemporaryMatrix<T>
1168 operator+(
const TemporaryMatrix<T> &a, T b)
1170 return const_cast<TemporaryMatrix<T> &
>(a) += b;
1180 template <
class T,
class C>
1181 inline TemporaryMatrix<T>
1184 return TemporaryMatrix<T>(b) += a;
1188 inline TemporaryMatrix<T>
1189 operator+(T a,
const TemporaryMatrix<T> &b)
1191 return const_cast<TemporaryMatrix<T> &
>(b) += a;
1201 template <
class T,
class C1,
class C2,
class C3>
1209 "subtract(): Matrix shapes must agree.");
1213 r(j, i) = a(j, i) - b(j, i);
1226 template <
class T,
class C1,
class C2>
1227 inline TemporaryMatrix<T>
1230 return TemporaryMatrix<T>(a) -= b;
1233 template <
class T,
class C>
1234 inline TemporaryMatrix<T>
1237 return const_cast<TemporaryMatrix<T> &
>(a) -= b;
1240 template <
class T,
class C>
1242 operator-(
const MultiArrayView<2, T, C> &a,
const TemporaryMatrix<T> &b)
1246 vigra_precondition(rows == b.rowCount() && cols == b.columnCount(),
1247 "Matrix::operator-(): Shape mismatch.");
1251 const_cast<TemporaryMatrix<T> &
>(b)(j, i) = a(j, i) - b(j, i);
1256 inline TemporaryMatrix<T>
1257 operator-(
const TemporaryMatrix<T> &a,
const TemporaryMatrix<T> &b)
1259 return const_cast<TemporaryMatrix<T> &
>(a) -= b;
1269 template <
class T,
class C>
1270 inline TemporaryMatrix<T>
1273 return TemporaryMatrix<T>(a) *= -NumericTraits<T>::one();
1277 inline TemporaryMatrix<T>
1280 return const_cast<TemporaryMatrix<T> &
>(a) *= -NumericTraits<T>::one();
1290 template <
class T,
class C>
1291 inline TemporaryMatrix<T>
1294 return TemporaryMatrix<T>(a) -= b;
1298 inline TemporaryMatrix<T>
1299 operator-(
const TemporaryMatrix<T> &a, T b)
1301 return const_cast<TemporaryMatrix<T> &
>(a) -= b;
1311 template <
class T,
class C>
1312 inline TemporaryMatrix<T>
1315 return TemporaryMatrix<T>(b.
shape(), a) -= b;
1330 template <
class T,
class C1,
class C2>
1331 typename NormTraits<T>::SquaredNormType
1334 typename NormTraits<T>::SquaredNormType ret =
1335 NumericTraits<typename NormTraits<T>::SquaredNormType>::zero();
1338 std::ptrdiff_t size = y.
shape(0);
1340 for(std::ptrdiff_t i = 0; i < size; ++i)
1341 ret += x(0, i) * y(i, 0);
1342 else if(x.
shape(1) == 1u && x.
shape(0) == size)
1343 for(std::ptrdiff_t i = 0; i < size; ++i)
1344 ret += x(i, 0) * y(i, 0);
1346 vigra_precondition(
false,
"dot(): wrong matrix shapes.");
1348 else if(y.
shape(0) == 1)
1350 std::ptrdiff_t size = y.
shape(1);
1352 for(std::ptrdiff_t i = 0; i < size; ++i)
1353 ret += x(0, i) * y(0, i);
1354 else if(x.
shape(1) == 1u && x.
shape(0) == size)
1355 for(std::ptrdiff_t i = 0; i < size; ++i)
1356 ret += x(i, 0) * y(0, i);
1358 vigra_precondition(
false,
"dot(): wrong matrix shapes.");
1361 vigra_precondition(
false,
"dot(): wrong matrix shapes.");
1372 template <
class T,
class C1,
class C2>
1373 typename NormTraits<T>::SquaredNormType
1378 "dot(): shape mismatch.");
1379 typename NormTraits<T>::SquaredNormType ret =
1380 NumericTraits<typename NormTraits<T>::SquaredNormType>::zero();
1393 template <
class T,
class C1,
class C2,
class C3>
1398 "cross(): vectors must have length 3.");
1399 r(0) = x(1)*y(2) - x(2)*y(1);
1400 r(1) = x(2)*y(0) - x(0)*y(2);
1401 r(2) = x(0)*y(1) - x(1)*y(0);
1412 template <
class T,
class C1,
class C2,
class C3>
1417 "cross(): vectors must have length 3.");
1418 r(0,0) = x(1,0)*y(2,0) - x(2,0)*y(1,0);
1419 r(1,0) = x(2,0)*y(0,0) - x(0,0)*y(2,0);
1420 r(2,0) = x(0,0)*y(1,0) - x(1,0)*y(0,0);
1431 template <
class T,
class C1,
class C2>
1435 TemporaryMatrix<T> ret(3, 1);
1448 template <
class T,
class C1,
class C2,
class C3>
1456 "outer(): shape mismatch.");
1459 r(j, i) = x(j, 0) * y(0, i);
1471 template <
class T,
class C1,
class C2>
1478 "outer(): shape mismatch.");
1479 TemporaryMatrix<T> ret(rows, cols);
1491 template <
class T,
class C>
1497 vigra_precondition(rows == 1 || cols == 1,
1498 "outer(): matrix does not represent a vector.");
1500 TemporaryMatrix<T> ret(size, size);
1506 ret(j, i) = x(0, j) * x(0, i);
1512 ret(j, i) = x(j, 0) * x(i, 0);
1523 PointWise(T
const & it)
1529 PointWise<T> pointWise(T
const & t)
1531 return PointWise<T>(t);
1542 template <
class T,
class C1,
class C2>
1548 "smul(): Matrix sizes must agree.");
1552 r(j, i) = a(j, i) * b;
1562 template <
class T,
class C2,
class C3>
1575 template <
class T,
class C1,
class C2,
class C3>
1583 "mmul(): Matrix shapes must agree.");
1589 r(j, i) = a(j, 0) * b(0, i);
1592 r(j, i) += a(j, k) * b(k, i);
1604 template <
class T,
class C1,
class C2>
1605 inline TemporaryMatrix<T>
1620 template <
class T,
class C1,
class C2,
class C3>
1628 "pmul(): Matrix shapes must agree.");
1632 r(j, i) = a(j, i) * b(j, i);
1645 template <
class T,
class C1,
class C2>
1646 inline TemporaryMatrix<T>
1649 TemporaryMatrix<T> ret(a.
shape());
1672 template <
class T,
class C,
class U>
1673 inline TemporaryMatrix<T>
1676 return pmul(a, b.t);
1686 template <
class T,
class C>
1687 inline TemporaryMatrix<T>
1690 return TemporaryMatrix<T>(a) *= b;
1694 inline TemporaryMatrix<T>
1695 operator*(
const TemporaryMatrix<T> &a, T b)
1697 return const_cast<TemporaryMatrix<T> &
>(a) *= b;
1707 template <
class T,
class C>
1708 inline TemporaryMatrix<T>
1711 return TemporaryMatrix<T>(b) *= a;
1715 inline TemporaryMatrix<T>
1716 operator*(T a,
const TemporaryMatrix<T> &b)
1718 return const_cast<TemporaryMatrix<T> &
>(b) *= a;
1729 template <
class T,
class A,
int N,
class DATA,
class DERIVED>
1734 "operator*(Matrix, TinyVector): Shape mismatch.");
1750 template <
class T,
int N,
class DATA,
class DERIVED,
class A>
1755 "operator*(TinyVector, Matrix): Shape mismatch.");
1771 template <
class T,
class C1,
class C2>
1772 inline TemporaryMatrix<T>
1787 template <
class T,
class C1,
class C2>
1793 "sdiv(): Matrix sizes must agree.");
1797 r(j, i) = a(j, i) / b;
1807 template <
class T,
class C1,
class C2,
class C3>
1815 "pdiv(): Matrix shapes must agree.");
1819 r(j, i) = a(j, i) / b(j, i);
1832 template <
class T,
class C1,
class C2>
1833 inline TemporaryMatrix<T>
1836 TemporaryMatrix<T> ret(a.
shape());
1859 template <
class T,
class C,
class U>
1860 inline TemporaryMatrix<T>
1863 return pdiv(a, b.t);
1873 template <
class T,
class C>
1874 inline TemporaryMatrix<T>
1877 return TemporaryMatrix<T>(a) /= b;
1881 inline TemporaryMatrix<T>
1882 operator/(
const TemporaryMatrix<T> &a, T b)
1884 return const_cast<TemporaryMatrix<T> &
>(a) /= b;
1894 template <
class T,
class C>
1895 inline TemporaryMatrix<T>
1898 return TemporaryMatrix<T>(b.
shape(), a) / pointWise(b);
1923 template <
class T,
class C>
1928 for(
int k=0; k < a.
size(); ++k)
1956 template <
class T,
class C>
1961 for(
int k=0; k < a.
size(); ++k)
1991 template <
class T,
class C,
class UnaryFunctor>
1996 for(
int k=0; k < a.
size(); ++k)
1998 if(condition(a[k]) && a[k] < vopt)
2026 template <
class T,
class C,
class UnaryFunctor>
2031 for(
int k=0; k < a.
size(); ++k)
2033 if(condition(a[k]) && vopt < a[k])
2044 template <
class T,
class C>
2047 linalg::TemporaryMatrix<T> t(v.
shape());
2057 linalg::TemporaryMatrix<T>
pow(linalg::TemporaryMatrix<T>
const & v, T exponent)
2059 linalg::TemporaryMatrix<T> & t =
const_cast<linalg::TemporaryMatrix<T> &
>(v);
2068 template <
class T,
class C>
2069 linalg::TemporaryMatrix<T>
pow(MultiArrayView<2, T, C>
const & v,
int exponent)
2071 linalg::TemporaryMatrix<T> t(v.shape());
2081 linalg::TemporaryMatrix<T>
pow(linalg::TemporaryMatrix<T>
const & v,
int exponent)
2083 linalg::TemporaryMatrix<T> & t =
const_cast<linalg::TemporaryMatrix<T> &
>(v);
2093 linalg::TemporaryMatrix<int>
pow(MultiArrayView<2, int, C>
const & v,
int exponent)
2095 linalg::TemporaryMatrix<int> t(v.shape());
2100 t(j, i) = (int)vigra::pow((
double)v(j, i), exponent);
2105 linalg::TemporaryMatrix<int>
pow(linalg::TemporaryMatrix<int>
const & v,
int exponent)
2107 linalg::TemporaryMatrix<int> & t =
const_cast<linalg::TemporaryMatrix<int> &
>(v);
2112 t(j, i) = (int)vigra::pow((
double)t(j, i), exponent);
2117 template <
class T,
class C>
2118 linalg::TemporaryMatrix<T>
sqrt(MultiArrayView<2, T, C>
const & v);
2120 template <
class T,
class C>
2121 linalg::TemporaryMatrix<T>
exp(MultiArrayView<2, T, C>
const & v);
2123 template <
class T,
class C>
2124 linalg::TemporaryMatrix<T>
log(MultiArrayView<2, T, C>
const & v);
2126 template <
class T,
class C>
2127 linalg::TemporaryMatrix<T>
log10(MultiArrayView<2, T, C>
const & v);
2129 template <
class T,
class C>
2130 linalg::TemporaryMatrix<T>
sin(MultiArrayView<2, T, C>
const & v);
2132 template <
class T,
class C>
2133 linalg::TemporaryMatrix<T>
asin(MultiArrayView<2, T, C>
const & v);
2135 template <
class T,
class C>
2136 linalg::TemporaryMatrix<T>
cos(MultiArrayView<2, T, C>
const & v);
2138 template <
class T,
class C>
2139 linalg::TemporaryMatrix<T>
acos(MultiArrayView<2, T, C>
const & v);
2141 template <
class T,
class C>
2142 linalg::TemporaryMatrix<T>
tan(MultiArrayView<2, T, C>
const & v);
2144 template <
class T,
class C>
2145 linalg::TemporaryMatrix<T>
atan(MultiArrayView<2, T, C>
const & v);
2147 template <
class T,
class C>
2148 linalg::TemporaryMatrix<T>
round(MultiArrayView<2, T, C>
const & v);
2150 template <
class T,
class C>
2151 linalg::TemporaryMatrix<T>
floor(MultiArrayView<2, T, C>
const & v);
2153 template <
class T,
class C>
2154 linalg::TemporaryMatrix<T>
ceil(MultiArrayView<2, T, C>
const & v);
2156 template <
class T,
class C>
2157 linalg::TemporaryMatrix<T>
abs(MultiArrayView<2, T, C>
const & v);
2159 template <
class T,
class C>
2160 linalg::TemporaryMatrix<T>
sq(MultiArrayView<2, T, C>
const & v);
2162 template <
class T,
class C>
2163 linalg::TemporaryMatrix<T>
sign(MultiArrayView<2, T, C>
const & v);
2165 #define VIGRA_MATRIX_UNARY_FUNCTION(FUNCTION, NAMESPACE) \
2166 using NAMESPACE::FUNCTION; \
2167 template <class T, class C> \
2168 linalg::TemporaryMatrix<T> FUNCTION(MultiArrayView<2, T, C> const & v) \
2170 linalg::TemporaryMatrix<T> t(v.shape()); \
2171 MultiArrayIndex m = rowCount(v), n = columnCount(v); \
2173 for(MultiArrayIndex i = 0; i < n; ++i) \
2174 for(MultiArrayIndex j = 0; j < m; ++j) \
2175 t(j, i) = NAMESPACE::FUNCTION(v(j, i)); \
2179 template <class T> \
2180 linalg::TemporaryMatrix<T> FUNCTION(linalg::Matrix<T> const & v) \
2182 linalg::TemporaryMatrix<T> t(v.shape()); \
2183 MultiArrayIndex m = rowCount(v), n = columnCount(v); \
2185 for(MultiArrayIndex i = 0; i < n; ++i) \
2186 for(MultiArrayIndex j = 0; j < m; ++j) \
2187 t(j, i) = NAMESPACE::FUNCTION(v(j, i)); \
2191 template <class T> \
2192 linalg::TemporaryMatrix<T> FUNCTION(linalg::TemporaryMatrix<T> const & v) \
2194 linalg::TemporaryMatrix<T> & t = const_cast<linalg::TemporaryMatrix<T> &>(v); \
2195 MultiArrayIndex m = rowCount(t), n = columnCount(t); \
2197 for(MultiArrayIndex i = 0; i < n; ++i) \
2198 for(MultiArrayIndex j = 0; j < m; ++j) \
2199 t(j, i) = NAMESPACE::FUNCTION(t(j, i)); \
2203 using linalg::FUNCTION;\
2206 VIGRA_MATRIX_UNARY_FUNCTION(
sqrt, std)
2207 VIGRA_MATRIX_UNARY_FUNCTION(
exp, std)
2208 VIGRA_MATRIX_UNARY_FUNCTION(
log, std)
2209 VIGRA_MATRIX_UNARY_FUNCTION(
log10, std)
2210 VIGRA_MATRIX_UNARY_FUNCTION(
sin, std)
2211 VIGRA_MATRIX_UNARY_FUNCTION(
asin, std)
2212 VIGRA_MATRIX_UNARY_FUNCTION(
cos, std)
2213 VIGRA_MATRIX_UNARY_FUNCTION(
acos, std)
2214 VIGRA_MATRIX_UNARY_FUNCTION(
tan, std)
2215 VIGRA_MATRIX_UNARY_FUNCTION(
atan, std)
2216 VIGRA_MATRIX_UNARY_FUNCTION(
round, vigra)
2217 VIGRA_MATRIX_UNARY_FUNCTION(
floor, vigra)
2218 VIGRA_MATRIX_UNARY_FUNCTION(
ceil, vigra)
2219 VIGRA_MATRIX_UNARY_FUNCTION(
abs, vigra)
2220 VIGRA_MATRIX_UNARY_FUNCTION(
sq, vigra)
2221 VIGRA_MATRIX_UNARY_FUNCTION(
sign, vigra)
2223 #undef VIGRA_MATRIX_UNARY_FUNCTION
2229 using linalg::RowMajor;
2230 using linalg::ColumnMajor;
2231 using linalg::Matrix;
2235 using linalg::pointWise;
2259 template <
class T,
class ALLOC>
2260 struct NormTraits<Matrix<T, ALLOC> >
2261 :
public NormTraits<MultiArray<2, T, ALLOC> >
2263 typedef NormTraits<MultiArray<2, T, ALLOC> > BaseType;
2264 typedef Matrix<T, ALLOC> Type;
2265 typedef typename BaseType::SquaredNormType SquaredNormType;
2266 typedef typename BaseType::NormType NormType;
2269 template <
class T,
class ALLOC>
2270 struct NormTraits<linalg::TemporaryMatrix<T, ALLOC> >
2271 :
public NormTraits<Matrix<T, ALLOC> >
2273 typedef NormTraits<Matrix<T, ALLOC> > BaseType;
2274 typedef linalg::TemporaryMatrix<T, ALLOC> Type;
2275 typedef typename BaseType::SquaredNormType SquaredNormType;
2276 typedef typename BaseType::NormType NormType;
2293 template <
class T,
class C>
2295 operator<<(ostream & s, const vigra::MultiArrayView<2, T, C> &m)
2299 ios::fmtflags flags = s.setf(ios::right | ios::fixed, ios::adjustfield | ios::floatfield);
2304 s << m(j, i) <<
" ";
2322 template <
class T1,
class C1,
class T2,
class C2,
class T3,
class C3>
2324 columnStatisticsImpl(MultiArrayView<2, T1, C1>
const & A,
2325 MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & sumOfSquaredDifferences)
2331 "columnStatistics(): Shape mismatch between input and output.");
2334 mean.init(NumericTraits<T2>::zero());
2335 sumOfSquaredDifferences.init(NumericTraits<T3>::zero());
2339 typedef typename NumericTraits<T2>::RealPromote TmpType;
2341 TmpType f = TmpType(1.0 / (k + 1.0)),
2342 f1 = TmpType(1.0 - f);
2344 sumOfSquaredDifferences += f1*
sq(t);
2348 template <
class T1,
class C1,
class T2,
class C2,
class T3,
class C3>
2350 columnStatistics2PassImpl(MultiArrayView<2, T1, C1>
const & A,
2351 MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & sumOfSquaredDifferences)
2357 "columnStatistics(): Shape mismatch between input and output.");
2360 mean.init(NumericTraits<T2>::zero());
2367 sumOfSquaredDifferences.init(NumericTraits<T3>::zero());
2370 sumOfSquaredDifferences +=
sq(
rowVector(A, k) - mean);
2435 template <
class T1,
class C1,
class T2,
class C2>
2438 MultiArrayView<2, T2, C2> & mean)
2443 "columnStatistics(): Shape mismatch between input and output.");
2445 mean.init(NumericTraits<T2>::zero());
2454 template <
class T1,
class C1,
class T2,
class C2,
class T3,
class C3>
2457 MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & stdDev)
2459 detail::columnStatisticsImpl(A, mean, stdDev);
2465 template <
class T1,
class C1,
class T2,
class C2,
class T3,
class C3,
class T4,
class C4>
2468 MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & stdDev, MultiArrayView<2, T4, C4> &
norm)
2475 "columnStatistics(): Shape mismatch between input and output.");
2477 detail::columnStatisticsImpl(A, mean, stdDev);
2478 norm =
sqrt(stdDev + T2(m) *
sq(mean));
2479 stdDev =
sqrt(stdDev / T3(m - 1.0));
2536 doxygen_overloaded_function(template <...>
void rowStatistics)
2538 template <
class T1,
class C1,
class T2,
class C2>
2541 MultiArrayView<2, T2, C2> & mean)
2544 "rowStatistics(): Shape mismatch between input and output.");
2545 MultiArrayView<2, T2, StridedArrayTag> tm =
transpose(mean);
2549 template <
class T1,
class C1,
class T2,
class C2,
class T3,
class C3>
2552 MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & stdDev)
2556 "rowStatistics(): Shape mismatch between input and output.");
2557 MultiArrayView<2, T2, StridedArrayTag> tm =
transpose(mean);
2558 MultiArrayView<2, T3, StridedArrayTag> ts =
transpose(stdDev);
2562 template <
class T1,
class C1,
class T2,
class C2,
class T3,
class C3,
class T4,
class C4>
2565 MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & stdDev, MultiArrayView<2, T4, C4> & norm)
2570 "rowStatistics(): Shape mismatch between input and output.");
2571 MultiArrayView<2, T2, StridedArrayTag> tm =
transpose(mean);
2572 MultiArrayView<2, T3, StridedArrayTag> ts =
transpose(stdDev);
2573 MultiArrayView<2, T4, StridedArrayTag> tn =
transpose(norm);
2579 template <
class T1,
class C1,
class U,
class T2,
class C2,
class T3,
class C3>
2580 void updateCovarianceMatrix(MultiArrayView<2, T1, C1>
const & features,
2581 U & count, MultiArrayView<2, T2, C2> & mean, MultiArrayView<2, T3, C3> & covariance)
2585 "updateCovarianceMatrix(): Features must be a row or column vector.");
2586 vigra_precondition(mean.shape() == features.shape(),
2587 "updateCovarianceMatrix(): Shape mismatch between feature vector and mean vector.");
2589 "updateCovarianceMatrix(): Shape mismatch between feature vector and covariance matrix.");
2592 Matrix<T2> t = features - mean;
2594 double f = 1.0 / count,
2602 covariance(k, k) += f1*
sq(t(0, k));
2605 covariance(k, l) += f1*t(0, k)*t(0, l);
2606 covariance(l, k) = covariance(k, l);
2614 covariance(k, k) += f1*
sq(t(k, 0));
2617 covariance(k, l) += f1*t(k, 0)*t(l, 0);
2618 covariance(l, k) = covariance(k, l);
2634 template <
class T1,
class C1,
class T2,
class C2>
2640 "covarianceMatrixOfColumns(): Shape mismatch between feature matrix and covariance matrix.");
2643 covariance.
init(NumericTraits<T2>::zero());
2645 detail::updateCovarianceMatrix(
rowVector(features, k), count, means, covariance);
2646 covariance /= T2(m - 1);
2657 template <
class T,
class C>
2674 template <
class T1,
class C1,
class T2,
class C2>
2680 "covarianceMatrixOfRows(): Shape mismatch between feature matrix and covariance matrix.");
2683 covariance.
init(NumericTraits<T2>::zero());
2685 detail::updateCovarianceMatrix(
columnVector(features, k), count, means, covariance);
2686 covariance /= T2(m - 1);
2697 template <
class T,
class C>
2706 enum DataPreparationGoals { ZeroMean = 1, UnitVariance = 2, UnitNorm = 4, UnitSum = 8 };
2708 inline DataPreparationGoals operator|(DataPreparationGoals l, DataPreparationGoals r)
2710 return DataPreparationGoals(
int(l) |
int(r));
2715 template <
class T,
class C1,
class C2,
class C3,
class C4>
2717 prepareDataImpl(
const MultiArrayView<2, T, C1> & A,
2718 MultiArrayView<2, T, C2> & res, MultiArrayView<2, T, C3> & offset, MultiArrayView<2, T, C4> & scaling,
2719 DataPreparationGoals goals)
2723 vigra_precondition(A.shape() == res.shape() &&
2726 "prepareDataImpl(): Shape mismatch between input and output.");
2731 offset.init(NumericTraits<T>::zero());
2732 scaling.init(NumericTraits<T>::one());
2736 bool zeroMean = (goals & ZeroMean) != 0;
2737 bool unitVariance = (goals & UnitVariance) != 0;
2738 bool unitNorm = (goals & UnitNorm) != 0;
2739 bool unitSum = (goals & UnitSum) != 0;
2743 vigra_precondition(goals == UnitSum,
2744 "prepareData(): Unit sum is not compatible with any other data preparation goal.");
2748 offset.init(NumericTraits<T>::zero());
2752 if(scaling(0, k) != NumericTraits<T>::zero())
2754 scaling(0, k) = NumericTraits<T>::one() / scaling(0, k);
2759 scaling(0, k) = NumericTraits<T>::one();
2766 vigra_precondition(!(unitVariance && unitNorm),
2767 "prepareData(): Unit variance and unit norm cannot be achieved at the same time.");
2769 Matrix<T> mean(1, n), sumOfSquaredDifferences(1, n);
2770 detail::columnStatisticsImpl(A, mean, sumOfSquaredDifferences);
2774 T stdDev =
std::sqrt(sumOfSquaredDifferences(0, k) / T(m-1));
2776 stdDev = NumericTraits<T>::zero();
2777 if(zeroMean && stdDev > NumericTraits<T>::zero())
2780 offset(0, k) = mean(0, k);
2781 mean(0, k) = NumericTraits<T>::zero();
2786 offset(0, k) = NumericTraits<T>::zero();
2789 T norm = mean(0,k) == NumericTraits<T>::zero()
2790 ?
std::sqrt(sumOfSquaredDifferences(0, k))
2791 : std::
sqrt(sumOfSquaredDifferences(0, k) + T(m) *
sq(mean(0,k)));
2792 if(unitNorm && norm > NumericTraits<T>::zero())
2795 scaling(0, k) = NumericTraits<T>::one() /
norm;
2797 else if(unitVariance && stdDev > NumericTraits<T>::zero())
2800 scaling(0, k) = NumericTraits<T>::one() / stdDev;
2804 scaling(0, k) = NumericTraits<T>::one();
2886 template <
class T,
class C1,
class C2,
class C3,
class C4>
2889 MultiArrayView<2, T, C2> & res, MultiArrayView<2, T, C3> & offset, MultiArrayView<2, T, C4> & scaling,
2890 DataPreparationGoals goals = ZeroMean | UnitVariance)
2892 detail::prepareDataImpl(A, res, offset, scaling, goals);
2895 template <
class T,
class C1,
class C2>
2897 prepareColumns(MultiArrayView<2, T, C1>
const & A, MultiArrayView<2, T, C2> & res,
2898 DataPreparationGoals goals = ZeroMean | UnitVariance)
2901 detail::prepareDataImpl(A, res, offset, scaling, goals);
2960 doxygen_overloaded_function(template <...>
void prepareRows)
2962 template <
class T,
class C1,
class C2,
class C3,
class C4>
2965 MultiArrayView<2, T, C2> & res, MultiArrayView<2, T, C3> & offset, MultiArrayView<2, T, C4> & scaling,
2966 DataPreparationGoals goals = ZeroMean | UnitVariance)
2969 detail::prepareDataImpl(
transpose(A), tr, to, ts, goals);
2972 template <
class T,
class C1,
class C2>
2974 prepareRows(MultiArrayView<2, T, C1>
const & A, MultiArrayView<2, T, C2> & res,
2975 DataPreparationGoals goals = ZeroMean | UnitVariance)
2977 MultiArrayView<2, T, StridedArrayTag> tr =
transpose(res);
2979 detail::prepareDataImpl(
transpose(A), tr, offset, scaling, goals);
2990 using linalg::ZeroMean;
2991 using linalg::UnitVariance;
2992 using linalg::UnitNorm;
2993 using linalg::UnitSum;
2999 #endif // VIGRA_MATRIX_HXX