core: cv::eigenNonSymmetric() via EigenvalueDecomposition

modules/core/include/opencv2/core.hpp

@@ -1913,8 +1913,9 @@ matrix src:
 @code
     src*eigenvectors.row(i).t() = eigenvalues.at<srcType>(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<srcType>(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

modules/core/src/lda.cpp

@@ -863,45 +863,49 @@ private:
         d = alloc_1d<double> (n);
         e = alloc_1d<double> (n);
         ort = alloc_1d<double> (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<double> (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<double> (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<double> (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<double> (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<double> eigenvalues64f;
+    eigensystem.eigenvalues().copyTo(eigenvalues64f);
+    CV_Assert(eigenvalues64f.size() == n);
+
+    std::vector<int> sort_indexes(n);
+    cv::sortIdx(eigenvalues64f, sort_indexes, SORT_EVERY_ROW | SORT_DESCENDING);
+
+    std::vector<double> 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_<double> sorted_eigenvectors64f((int)n, (int)n, CV_64FC1);
+        for (size_t i = 0; i < n; i++)
+        {
+            double* pDst = sorted_eigenvectors64f.ptr<double>((int)i);
+            double* pSrc = eigenvectors64f.ptr<double>(sort_indexes[(int)i]);
+            CV_Assert(pSrc != NULL);
+            memcpy(pDst, pSrc, n * sizeof(double));
+        }
+        sorted_eigenvectors64f.convertTo(_evects, type);
+    }
+}
+
 
 //------------------------------------------------------------------------------
 // Linear Discriminant Analysis implementation

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<typename T>
+static void testEigen(const Mat_<T>& src, const Mat_<T>& expected_eigenvalues, bool runSymmetric = false)
+{
+    SCOPED_TRACE(runSymmetric ? "cv::eigen" : "cv::eigenNonSymmetric");
+
+    int type = traits::Type<T>::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<srcType>(i)*eigenvectors.row(i).t()
+    for (int i = 0; i < src.rows; i++)
+    {
+        EXPECT_NEAR(eigenvalues.at<T>(i), expected_eigenvalues(i), eps) << "i=" << i;
+        Mat lhs = src*eigenvectors.row(i).t();
+        Mat rhs = eigenvalues.at<T>(i)*eigenvectors.row(i).t();
+        EXPECT_LT(cvtest::norm(lhs, rhs, NORM_INF), eps)
+                << "i=" << i << " eigenvalue=" << eigenvalues.at<T>(i) << std::endl
+                << "lhs=" << lhs.t() << std::endl
+                << "rhs=" << rhs.t();
+    }
+}
+
+template<typename T>
+static void testEigenSymmetric3x3()
+{
+    /*const*/ T values_[] = {
+            2, -1, 0,
+            -1, 2, -1,
+            0, -1, 2
+    };
+    Mat_<T> src(3, 3, values_);
+
+    /*const*/ T expected_eigenvalues_[] = { 3.414213562373095f, 2, 0.585786437626905f };
+    Mat_<T> expected_eigenvalues(3, 1, expected_eigenvalues_);
+
+    testEigen(src, expected_eigenvalues);
+    testEigen(src, expected_eigenvalues, true);
+}
+TEST(Core_EigenSymmetric, float3x3) { testEigenSymmetric3x3<float>(); }
+TEST(Core_EigenSymmetric, double3x3) { testEigenSymmetric3x3<double>(); }
+
+template<typename T>
+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_<T> src(5, 5, values_);
+
+    /*const*/ T expected_eigenvalues_[] = { 7.028919644935684f, 4.406130784616501f, 3.73626552682258f, 1.438067799899037f, 0.390616243726198f };
+    Mat_<T> expected_eigenvalues(5, 1, expected_eigenvalues_);
+
+    testEigen(src, expected_eigenvalues);
+    testEigen(src, expected_eigenvalues, true);
+}
+TEST(Core_EigenSymmetric, float5x5) { testEigenSymmetric5x5<float>(); }
+TEST(Core_EigenSymmetric, double5x5) { testEigenSymmetric5x5<double>(); }
+
+
+template<typename T>
+static void testEigen2x2()
+{
+    /*const*/ T values_[] = { 4, 1, 6, 3 };
+    Mat_<T> src(2, 2, values_);
+
+    /*const*/ T expected_eigenvalues_[] = { 6, 1 };
+    Mat_<T> expected_eigenvalues(2, 1, expected_eigenvalues_);
+
+    testEigen(src, expected_eigenvalues);
+}
+TEST(Core_EigenNonSymmetric, float2x2) { testEigen2x2<float>(); }
+TEST(Core_EigenNonSymmetric, double2x2) { testEigen2x2<double>(); }
+
+template<typename T>
+static void testEigen3x3()
+{
+    /*const*/ T values_[3*3] = {
+            3,1,0,
+            0,3,1,
+            0,0,3
+    };
+    Mat_<T> src(3, 3, values_);
+
+    /*const*/ T expected_eigenvalues_[] = { 3, 3, 3 };
+    Mat_<T> expected_eigenvalues(3, 1, expected_eigenvalues_);
+
+    testEigen(src, expected_eigenvalues);
+}
+TEST(Core_EigenNonSymmetric, float3x3) { testEigen3x3<float>(); }
+TEST(Core_EigenNonSymmetric, double3x3) { testEigen3x3<double>(); }