From 529632f8d000c36e9427010dbe73d84eaefd8c20 Mon Sep 17 00:00:00 2001 From: Alexander Alekhin Date: Sat, 30 Sep 2017 11:55:13 +0000 Subject: [PATCH] core: cv::eigenNonSymmetric() via EigenvalueDecomposition --- modules/core/include/opencv2/core.hpp | 24 ++++- modules/core/src/lda.cpp | 101 ++++++++++++++++----- modules/core/test/test_eigen.cpp | 121 ++++++++++++++++++++++++++ 3 files changed, 219 insertions(+), 27 deletions(-) diff --git a/modules/core/include/opencv2/core.hpp b/modules/core/include/opencv2/core.hpp index 880c04c3e1..8054d31b28 100644 --- a/modules/core/include/opencv2/core.hpp +++ b/modules/core/include/opencv2/core.hpp @@ -1913,8 +1913,9 @@ matrix src: @code src*eigenvectors.row(i).t() = eigenvalues.at(i)*eigenvectors.row(i).t() @endcode -@note in the new and the old interfaces different ordering of eigenvalues and eigenvectors -parameters is used. + +@note Use cv::eigenNonSymmetric for calculation of real eigenvalues and eigenvectors of non-symmetric matrix. + @param src input matrix that must have CV_32FC1 or CV_64FC1 type, square size and be symmetrical (src ^T^ == src). @param eigenvalues output vector of eigenvalues of the same type as src; the eigenvalues are stored @@ -1922,11 +1923,28 @@ in the descending order. @param eigenvectors output matrix of eigenvectors; it has the same size and type as src; the eigenvectors are stored as subsequent matrix rows, in the same order as the corresponding eigenvalues. -@sa completeSymm , PCA +@sa eigenNonSymmetric, completeSymm , PCA */ CV_EXPORTS_W bool eigen(InputArray src, OutputArray eigenvalues, OutputArray eigenvectors = noArray()); +/** @brief Calculates eigenvalues and eigenvectors of a non-symmetric matrix (real eigenvalues only). + +@note Assumes real eigenvalues. + +The function calculates eigenvalues and eigenvectors (optional) of the square matrix src: +@code + src*eigenvectors.row(i).t() = eigenvalues.at(i)*eigenvectors.row(i).t() +@endcode + +@param src input matrix (CV_32FC1 or CV_64FC1 type). +@param eigenvalues output vector of eigenvalues (type is the same type as src). +@param eigenvectors output matrix of eigenvectors (type is the same type as src). The eigenvectors are stored as subsequent matrix rows, in the same order as the corresponding eigenvalues. +@sa eigen +*/ +CV_EXPORTS_W void eigenNonSymmetric(InputArray src, OutputArray eigenvalues, + OutputArray eigenvectors); + /** @brief Calculates the covariance matrix of a set of vectors. The function cv::calcCovarMatrix calculates the covariance matrix and, optionally, the mean vector of diff --git a/modules/core/src/lda.cpp b/modules/core/src/lda.cpp index 20f0604093..bf356ec53e 100644 --- a/modules/core/src/lda.cpp +++ b/modules/core/src/lda.cpp @@ -863,45 +863,49 @@ private: d = alloc_1d (n); e = alloc_1d (n); ort = alloc_1d (n); - // Reduce to Hessenberg form. - orthes(); - // Reduce Hessenberg to real Schur form. - hqr2(); - // Copy eigenvalues to OpenCV Matrix. - _eigenvalues.create(1, n, CV_64FC1); - for (int i = 0; i < n; i++) { - _eigenvalues.at (0, i) = d[i]; + try { + // Reduce to Hessenberg form. + orthes(); + // Reduce Hessenberg to real Schur form. + hqr2(); + // Copy eigenvalues to OpenCV Matrix. + _eigenvalues.create(1, n, CV_64FC1); + for (int i = 0; i < n; i++) { + _eigenvalues.at (0, i) = d[i]; + } + // Copy eigenvectors to OpenCV Matrix. + _eigenvectors.create(n, n, CV_64FC1); + for (int i = 0; i < n; i++) + for (int j = 0; j < n; j++) + _eigenvectors.at (i, j) = V[i][j]; + // Deallocate the memory by releasing all internal working data. + release(); + } + catch (...) + { + release(); + throw; } - // Copy eigenvectors to OpenCV Matrix. - _eigenvectors.create(n, n, CV_64FC1); - for (int i = 0; i < n; i++) - for (int j = 0; j < n; j++) - _eigenvectors.at (i, j) = V[i][j]; - // Deallocate the memory by releasing all internal working data. - release(); } public: - EigenvalueDecomposition() - : n(0), cdivr(0), cdivi(0), d(0), e(0), ort(0), V(0), H(0) {} - // Initializes & computes the Eigenvalue Decomposition for a general matrix // given in src. This function is a port of the EigenvalueSolver in JAMA, // which has been released to public domain by The MathWorks and the // National Institute of Standards and Technology (NIST). - EigenvalueDecomposition(InputArray src) { - compute(src); + EigenvalueDecomposition(InputArray src, bool fallbackSymmetric = true) { + compute(src, fallbackSymmetric); } // This function computes the Eigenvalue Decomposition for a general matrix // given in src. This function is a port of the EigenvalueSolver in JAMA, // which has been released to public domain by The MathWorks and the // National Institute of Standards and Technology (NIST). - void compute(InputArray src) + void compute(InputArray src, bool fallbackSymmetric) { CV_INSTRUMENT_REGION() - if(isSymmetric(src)) { + if(fallbackSymmetric && isSymmetric(src)) { // Fall back to OpenCV for a symmetric matrix! cv::eigen(src, _eigenvalues, _eigenvectors); } else { @@ -930,11 +934,60 @@ public: ~EigenvalueDecomposition() {} // Returns the eigenvalues of the Eigenvalue Decomposition. - Mat eigenvalues() { return _eigenvalues; } + Mat eigenvalues() const { return _eigenvalues; } // Returns the eigenvectors of the Eigenvalue Decomposition. - Mat eigenvectors() { return _eigenvectors; } + Mat eigenvectors() const { return _eigenvectors; } }; +void eigenNonSymmetric(InputArray _src, OutputArray _evals, OutputArray _evects) +{ + CV_INSTRUMENT_REGION() + + Mat src = _src.getMat(); + int type = src.type(); + size_t n = (size_t)src.rows; + + CV_Assert(src.rows == src.cols); + CV_Assert(type == CV_32F || type == CV_64F); + + Mat src64f; + if (type == CV_32F) + src.convertTo(src64f, CV_32FC1); + else + src64f = src; + + EigenvalueDecomposition eigensystem(src64f, false); + + // EigenvalueDecomposition returns transposed and non-sorted eigenvalues + std::vector eigenvalues64f; + eigensystem.eigenvalues().copyTo(eigenvalues64f); + CV_Assert(eigenvalues64f.size() == n); + + std::vector sort_indexes(n); + cv::sortIdx(eigenvalues64f, sort_indexes, SORT_EVERY_ROW | SORT_DESCENDING); + + std::vector sorted_eigenvalues64f(n); + for (size_t i = 0; i < n; i++) sorted_eigenvalues64f[i] = eigenvalues64f[sort_indexes[i]]; + + Mat(sorted_eigenvalues64f).convertTo(_evals, type); + + if( _evects.needed() ) + { + Mat eigenvectors64f = eigensystem.eigenvectors().t(); // transpose + CV_Assert((size_t)eigenvectors64f.rows == n); + CV_Assert((size_t)eigenvectors64f.cols == n); + Mat_ sorted_eigenvectors64f((int)n, (int)n, CV_64FC1); + for (size_t i = 0; i < n; i++) + { + double* pDst = sorted_eigenvectors64f.ptr((int)i); + double* pSrc = eigenvectors64f.ptr(sort_indexes[(int)i]); + CV_Assert(pSrc != NULL); + memcpy(pDst, pSrc, n * sizeof(double)); + } + sorted_eigenvectors64f.convertTo(_evects, type); + } +} + //------------------------------------------------------------------------------ // Linear Discriminant Analysis implementation diff --git a/modules/core/test/test_eigen.cpp b/modules/core/test/test_eigen.cpp index bd51c7441a..44f4a729c8 100644 --- a/modules/core/test/test_eigen.cpp +++ b/modules/core/test/test_eigen.cpp @@ -412,3 +412,124 @@ TEST(Core_Eigen, scalar_32) {Core_EigenTest_Scalar_32 test; test.safe_run(); } TEST(Core_Eigen, scalar_64) {Core_EigenTest_Scalar_64 test; test.safe_run(); } TEST(Core_Eigen, vector_32) { Core_EigenTest_32 test; test.safe_run(); } TEST(Core_Eigen, vector_64) { Core_EigenTest_64 test; test.safe_run(); } + +template +static void testEigen(const Mat_& src, const Mat_& expected_eigenvalues, bool runSymmetric = false) +{ + SCOPED_TRACE(runSymmetric ? "cv::eigen" : "cv::eigenNonSymmetric"); + + int type = traits::Type::value; + const T eps = 1e-6f; + + Mat eigenvalues, eigenvectors, eigenvalues0; + + if (runSymmetric) + { + cv::eigen(src, eigenvalues0, noArray()); + cv::eigen(src, eigenvalues, eigenvectors); + } + else + { + cv::eigenNonSymmetric(src, eigenvalues0, noArray()); + cv::eigenNonSymmetric(src, eigenvalues, eigenvectors); + } +#if 0 + std::cout << "src = " << src << std::endl; + std::cout << "eigenvalues.t() = " << eigenvalues.t() << std::endl; + std::cout << "eigenvectors = " << eigenvectors << std::endl; +#endif + ASSERT_EQ(type, eigenvalues0.type()); + ASSERT_EQ(type, eigenvalues.type()); + ASSERT_EQ(type, eigenvectors.type()); + + ASSERT_EQ(src.rows, eigenvalues.rows); + ASSERT_EQ(eigenvalues.rows, eigenvectors.rows); + ASSERT_EQ(src.rows, eigenvectors.cols); + + EXPECT_LT(cvtest::norm(eigenvalues, eigenvalues0, NORM_INF), eps); + + // check definition: src*eigenvectors.row(i).t() = eigenvalues.at(i)*eigenvectors.row(i).t() + for (int i = 0; i < src.rows; i++) + { + EXPECT_NEAR(eigenvalues.at(i), expected_eigenvalues(i), eps) << "i=" << i; + Mat lhs = src*eigenvectors.row(i).t(); + Mat rhs = eigenvalues.at(i)*eigenvectors.row(i).t(); + EXPECT_LT(cvtest::norm(lhs, rhs, NORM_INF), eps) + << "i=" << i << " eigenvalue=" << eigenvalues.at(i) << std::endl + << "lhs=" << lhs.t() << std::endl + << "rhs=" << rhs.t(); + } +} + +template +static void testEigenSymmetric3x3() +{ + /*const*/ T values_[] = { + 2, -1, 0, + -1, 2, -1, + 0, -1, 2 + }; + Mat_ src(3, 3, values_); + + /*const*/ T expected_eigenvalues_[] = { 3.414213562373095f, 2, 0.585786437626905f }; + Mat_ expected_eigenvalues(3, 1, expected_eigenvalues_); + + testEigen(src, expected_eigenvalues); + testEigen(src, expected_eigenvalues, true); +} +TEST(Core_EigenSymmetric, float3x3) { testEigenSymmetric3x3(); } +TEST(Core_EigenSymmetric, double3x3) { testEigenSymmetric3x3(); } + +template +static void testEigenSymmetric5x5() +{ + /*const*/ T values_[5*5] = { + 5, -1, 0, 2, 1, + -1, 4, -1, 0, 0, + 0, -1, 3, 1, -1, + 2, 0, 1, 4, 0, + 1, 0, -1, 0, 1 + }; + Mat_ src(5, 5, values_); + + /*const*/ T expected_eigenvalues_[] = { 7.028919644935684f, 4.406130784616501f, 3.73626552682258f, 1.438067799899037f, 0.390616243726198f }; + Mat_ expected_eigenvalues(5, 1, expected_eigenvalues_); + + testEigen(src, expected_eigenvalues); + testEigen(src, expected_eigenvalues, true); +} +TEST(Core_EigenSymmetric, float5x5) { testEigenSymmetric5x5(); } +TEST(Core_EigenSymmetric, double5x5) { testEigenSymmetric5x5(); } + + +template +static void testEigen2x2() +{ + /*const*/ T values_[] = { 4, 1, 6, 3 }; + Mat_ src(2, 2, values_); + + /*const*/ T expected_eigenvalues_[] = { 6, 1 }; + Mat_ expected_eigenvalues(2, 1, expected_eigenvalues_); + + testEigen(src, expected_eigenvalues); +} +TEST(Core_EigenNonSymmetric, float2x2) { testEigen2x2(); } +TEST(Core_EigenNonSymmetric, double2x2) { testEigen2x2(); } + +template +static void testEigen3x3() +{ + /*const*/ T values_[3*3] = { + 3,1,0, + 0,3,1, + 0,0,3 + }; + Mat_ src(3, 3, values_); + + /*const*/ T expected_eigenvalues_[] = { 3, 3, 3 }; + Mat_ expected_eigenvalues(3, 1, expected_eigenvalues_); + + testEigen(src, expected_eigenvalues); +} +TEST(Core_EigenNonSymmetric, float3x3) { testEigen3x3(); } +TEST(Core_EigenNonSymmetric, double3x3) { testEigen3x3(); }