@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
@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
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 {
~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
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>(); }
projects / platform / upstream / opencv.git / commitdiff