Use Eigen::SelfAdjointEigenSolver in cv::eigen

woody.chow · alalek · commit 611cf8d86f5a · 2017-12-05T02:40:55.000+03:00
diff --git a/modules/core/src/lapack.cpp b/modules/core/src/lapack.cpp
@@ -43,6 +43,12 @@
 #include "precomp.hpp"
 #include <limits>
 
+#ifdef HAVE_EIGEN
+#include <Eigen/Core>
+#include <Eigen/Eigenvalues>
+#include "opencv2/core/eigen.hpp"
+#endif
+
 #if defined _M_IX86 && defined _MSC_VER && _MSC_VER < 1700
 #pragma float_control(precise, on)
 #endif
@@ -1396,6 +1402,47 @@ bool cv::eigen( InputArray _src, OutputArray _evals, OutputArray _evects )
         v = _evects.getMat();
     }
 
+#ifdef HAVE_EIGEN
+    const bool evecNeeded = _evects.needed();
+    const int esOptions = evecNeeded ? Eigen::ComputeEigenvectors : Eigen::EigenvaluesOnly;
+    _evals.create(n, 1, type);
+    cv::Mat evals = _evals.getMat();
+    if ( type == CV_64F )
+    {
+        Eigen::MatrixXd src_eig, zeros_eig;
+        cv::cv2eigen(src, src_eig);
+
+        Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es;
+        es.compute(src_eig, esOptions);
+        if ( es.info() == Eigen::Success )
+        {
+            cv::eigen2cv(es.eigenvalues().reverse().eval(), evals);
+            if ( evecNeeded )
+            {
+                cv::Mat evects = _evects.getMat();
+                cv::eigen2cv(es.eigenvectors().rowwise().reverse().transpose().eval(), v);
+            }
+            return true;
+        }
+    } else { // CV_32F
+        Eigen::MatrixXf src_eig, zeros_eig;
+        cv::cv2eigen(src, src_eig);
+
+        Eigen::SelfAdjointEigenSolver<Eigen::MatrixXf> es;
+        es.compute(src_eig, esOptions);
+        if ( es.info() == Eigen::Success )
+        {
+            cv::eigen2cv(es.eigenvalues().reverse().eval(), evals);
+            if ( evecNeeded )
+            {
+                cv::eigen2cv(es.eigenvectors().rowwise().reverse().transpose().eval(), v);
+            }
+            return true;
+        }
+    }
+    return false;
+#else
+
     size_t elemSize = src.elemSize(), astep = alignSize(n*elemSize, 16);
     AutoBuffer<uchar> buf(n*astep + n*5*elemSize + 32);
     uchar* ptr = alignPtr((uchar*)buf, 16);
@@ -1408,6 +1455,7 @@ bool cv::eigen( InputArray _src, OutputArray _evals, OutputArray _evects )
 
     w.copyTo(_evals);
     return ok;
+#endif
 }
 
 namespace cv
diff --git a/modules/core/test/test_eigen.cpp b/modules/core/test/test_eigen.cpp
@@ -59,7 +59,7 @@ using namespace std;
 #define MESSAGE_ERROR_DIFF_1 "Accuracy of eigen values computing less than required."
 #define MESSAGE_ERROR_DIFF_2 "Accuracy of eigen vectors computing less than required."
 #define MESSAGE_ERROR_ORTHO "Matrix of eigen vectors is not orthogonal."
-#define MESSAGE_ERROR_ORDER "Eigen values are not sorted in ascending order."
+#define MESSAGE_ERROR_ORDER "Eigen values are not sorted in descending order."
 
 const int COUNT_NORM_TYPES = 3;
 const int NORM_TYPE[COUNT_NORM_TYPES] = {cv::NORM_L1, cv::NORM_L2, cv::NORM_INF};
@@ -257,7 +257,7 @@ bool Core_EigenTest::check_pairs_order(const cv::Mat& eigen_values)
                 if (!(eigen_values.at<float>(i, 0) > eigen_values.at<float>(i+1, 0)))
                 {
                 std::cout << endl; std::cout << "Checking order of eigen values vector " << eigen_values << "..." << endl;
-                std::cout << "Pair of indexes with non ascending of eigen values: (" << i << ", " << i+1 << ")." << endl;
+                std::cout << "Pair of indexes with non descending of eigen values: (" << i << ", " << i+1 << ")." << endl;
                 std::cout << endl;
                 CV_Error(CORE_EIGEN_ERROR_ORDER, MESSAGE_ERROR_ORDER);
                 return false;
@@ -272,9 +272,9 @@ bool Core_EigenTest::check_pairs_order(const cv::Mat& eigen_values)
                 if (!(eigen_values.at<double>(i, 0) > eigen_values.at<double>(i+1, 0)))
                 {
                     std::cout << endl; std::cout << "Checking order of eigen values vector " << eigen_values << "..." << endl;
-                    std::cout << "Pair of indexes with non ascending of eigen values: (" << i << ", " << i+1 << ")." << endl;
+                    std::cout << "Pair of indexes with non descending of eigen values: (" << i << ", " << i+1 << ")." << endl;
                     std::cout << endl;
-                    CV_Error(CORE_EIGEN_ERROR_ORDER, "Eigen values are not sorted in ascending order.");
+                    CV_Error(CORE_EIGEN_ERROR_ORDER, "Eigen values are not sorted in descending order.");
                     return false;
                 }
 
