Skip to content

Commit f646f61

Browse files
vpisarevvpisarev
committed
Merge pull request opencv#9752 from alalek:eigenNonSymmetric
2 parents ca9a3af + 529632f commit f646f61

File tree

3 files changed

+219
-27
lines changed

3 files changed

+219
-27
lines changed

‎modules/core/include/opencv2/core.hpp

Lines changed: 21 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1913,20 +1913,38 @@ matrix src:
19131913
@code
19141914
src*eigenvectors.row(i).t() = eigenvalues.at<srcType>(i)*eigenvectors.row(i).t()
19151915
@endcode
1916-
@note in the new and the old interfaces different ordering of eigenvalues and eigenvectors
1917-
parameters is used.
1916+
1917+
@note Use cv::eigenNonSymmetric for calculation of real eigenvalues and eigenvectors of non-symmetric matrix.
1918+
19181919
@param src input matrix that must have CV_32FC1 or CV_64FC1 type, square size and be symmetrical
19191920
(src ^T^ == src).
19201921
@param eigenvalues output vector of eigenvalues of the same type as src; the eigenvalues are stored
19211922
in the descending order.
19221923
@param eigenvectors output matrix of eigenvectors; it has the same size and type as src; the
19231924
eigenvectors are stored as subsequent matrix rows, in the same order as the corresponding
19241925
eigenvalues.
1925-
@sa completeSymm , PCA
1926+
@sa eigenNonSymmetric, completeSymm , PCA
19261927
*/
19271928
CV_EXPORTS_W bool eigen(InputArray src, OutputArray eigenvalues,
19281929
OutputArray eigenvectors = noArray());
19291930

1931+
/** @brief Calculates eigenvalues and eigenvectors of a non-symmetric matrix (real eigenvalues only).
1932+
1933+
@note Assumes real eigenvalues.
1934+
1935+
The function calculates eigenvalues and eigenvectors (optional) of the square matrix src:
1936+
@code
1937+
src*eigenvectors.row(i).t() = eigenvalues.at<srcType>(i)*eigenvectors.row(i).t()
1938+
@endcode
1939+
1940+
@param src input matrix (CV_32FC1 or CV_64FC1 type).
1941+
@param eigenvalues output vector of eigenvalues (type is the same type as src).
1942+
@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.
1943+
@sa eigen
1944+
*/
1945+
CV_EXPORTS_W void eigenNonSymmetric(InputArray src, OutputArray eigenvalues,
1946+
OutputArray eigenvectors);
1947+
19301948
/** @brief Calculates the covariance matrix of a set of vectors.
19311949
19321950
The function cv::calcCovarMatrix calculates the covariance matrix and, optionally, the mean vector of

‎modules/core/src/lda.cpp

Lines changed: 77 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -863,45 +863,49 @@ class EigenvalueDecomposition {
863863
d = alloc_1d<double> (n);
864864
e = alloc_1d<double> (n);
865865
ort = alloc_1d<double> (n);
866-
// Reduce to Hessenberg form.
867-
orthes();
868-
// Reduce Hessenberg to real Schur form.
869-
hqr2();
870-
// Copy eigenvalues to OpenCV Matrix.
871-
_eigenvalues.create(1, n, CV_64FC1);
872-
for (int i = 0; i < n; i++) {
873-
_eigenvalues.at<double> (0, i) = d[i];
866+
try {
867+
// Reduce to Hessenberg form.
868+
orthes();
869+
// Reduce Hessenberg to real Schur form.
870+
hqr2();
871+
// Copy eigenvalues to OpenCV Matrix.
872+
_eigenvalues.create(1, n, CV_64FC1);
873+
for (int i = 0; i < n; i++) {
874+
_eigenvalues.at<double> (0, i) = d[i];
875+
}
876+
// Copy eigenvectors to OpenCV Matrix.
877+
_eigenvectors.create(n, n, CV_64FC1);
878+
for (int i = 0; i < n; i++)
879+
for (int j = 0; j < n; j++)
880+
_eigenvectors.at<double> (i, j) = V[i][j];
881+
// Deallocate the memory by releasing all internal working data.
882+
release();
883+
}
884+
catch (...)
885+
{
886+
release();
887+
throw;
874888
}
875-
// Copy eigenvectors to OpenCV Matrix.
876-
_eigenvectors.create(n, n, CV_64FC1);
877-
for (int i = 0; i < n; i++)
878-
for (int j = 0; j < n; j++)
879-
_eigenvectors.at<double> (i, j) = V[i][j];
880-
// Deallocate the memory by releasing all internal working data.
881-
release();
882889
}
883890

884891
public:
885-
EigenvalueDecomposition()
886-
: n(0), cdivr(0), cdivi(0), d(0), e(0), ort(0), V(0), H(0) {}
887-
888892
// Initializes & computes the Eigenvalue Decomposition for a general matrix
889893
// given in src. This function is a port of the EigenvalueSolver in JAMA,
890894
// which has been released to public domain by The MathWorks and the
891895
// National Institute of Standards and Technology (NIST).
892-
EigenvalueDecomposition(InputArray src) {
893-
compute(src);
896+
EigenvalueDecomposition(InputArray src, bool fallbackSymmetric = true) {
897+
compute(src, fallbackSymmetric);
894898
}
895899

896900
// This function computes the Eigenvalue Decomposition for a general matrix
897901
// given in src. This function is a port of the EigenvalueSolver in JAMA,
898902
// which has been released to public domain by The MathWorks and the
899903
// National Institute of Standards and Technology (NIST).
900-
void compute(InputArray src)
904+
void compute(InputArray src, bool fallbackSymmetric)
901905
{
902906
CV_INSTRUMENT_REGION()
903907

904-
if(isSymmetric(src)) {
908+
if(fallbackSymmetric && isSymmetric(src)) {
905909
// Fall back to OpenCV for a symmetric matrix!
906910
cv::eigen(src, _eigenvalues, _eigenvectors);
907911
} else {
@@ -930,11 +934,60 @@ class EigenvalueDecomposition {
930934
~EigenvalueDecomposition() {}
931935

932936
// Returns the eigenvalues of the Eigenvalue Decomposition.
933-
Mat eigenvalues() { return _eigenvalues; }
937+
Mat eigenvalues() const { return _eigenvalues; }
934938
// Returns the eigenvectors of the Eigenvalue Decomposition.
935-
Mat eigenvectors() { return _eigenvectors; }
939+
Mat eigenvectors() const { return _eigenvectors; }
936940
};
937941

942+
void eigenNonSymmetric(InputArray _src, OutputArray _evals, OutputArray _evects)
943+
{
944+
CV_INSTRUMENT_REGION()
945+
946+
Mat src = _src.getMat();
947+
int type = src.type();
948+
size_t n = (size_t)src.rows;
949+
950+
CV_Assert(src.rows == src.cols);
951+
CV_Assert(type == CV_32F || type == CV_64F);
952+
953+
Mat src64f;
954+
if (type == CV_32F)
955+
src.convertTo(src64f, CV_32FC1);
956+
else
957+
src64f = src;
958+
959+
EigenvalueDecomposition eigensystem(src64f, false);
960+
961+
// EigenvalueDecomposition returns transposed and non-sorted eigenvalues
962+
std::vector<double> eigenvalues64f;
963+
eigensystem.eigenvalues().copyTo(eigenvalues64f);
964+
CV_Assert(eigenvalues64f.size() == n);
965+
966+
std::vector<int> sort_indexes(n);
967+
cv::sortIdx(eigenvalues64f, sort_indexes, SORT_EVERY_ROW | SORT_DESCENDING);
968+
969+
std::vector<double> sorted_eigenvalues64f(n);
970+
for (size_t i = 0; i < n; i++) sorted_eigenvalues64f[i] = eigenvalues64f[sort_indexes[i]];
971+
972+
Mat(sorted_eigenvalues64f).convertTo(_evals, type);
973+
974+
if( _evects.needed() )
975+
{
976+
Mat eigenvectors64f = eigensystem.eigenvectors().t(); // transpose
977+
CV_Assert((size_t)eigenvectors64f.rows == n);
978+
CV_Assert((size_t)eigenvectors64f.cols == n);
979+
Mat_<double> sorted_eigenvectors64f((int)n, (int)n, CV_64FC1);
980+
for (size_t i = 0; i < n; i++)
981+
{
982+
double* pDst = sorted_eigenvectors64f.ptr<double>((int)i);
983+
double* pSrc = eigenvectors64f.ptr<double>(sort_indexes[(int)i]);
984+
CV_Assert(pSrc != NULL);
985+
memcpy(pDst, pSrc, n * sizeof(double));
986+
}
987+
sorted_eigenvectors64f.convertTo(_evects, type);
988+
}
989+
}
990+
938991

939992
//------------------------------------------------------------------------------
940993
// Linear Discriminant Analysis implementation

‎modules/core/test/test_eigen.cpp

Lines changed: 121 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -412,3 +412,124 @@ TEST(Core_Eigen, scalar_32) {Core_EigenTest_Scalar_32 test; test.safe_run(); }
412412
TEST(Core_Eigen, scalar_64) {Core_EigenTest_Scalar_64 test; test.safe_run(); }
413413
TEST(Core_Eigen, vector_32) { Core_EigenTest_32 test; test.safe_run(); }
414414
TEST(Core_Eigen, vector_64) { Core_EigenTest_64 test; test.safe_run(); }
415+
416+
template<typename T>
417+
static void testEigen(const Mat_<T>& src, const Mat_<T>& expected_eigenvalues, bool runSymmetric = false)
418+
{
419+
SCOPED_TRACE(runSymmetric ? "cv::eigen" : "cv::eigenNonSymmetric");
420+
421+
int type = traits::Type<T>::value;
422+
const T eps = 1e-6f;
423+
424+
Mat eigenvalues, eigenvectors, eigenvalues0;
425+
426+
if (runSymmetric)
427+
{
428+
cv::eigen(src, eigenvalues0, noArray());
429+
cv::eigen(src, eigenvalues, eigenvectors);
430+
}
431+
else
432+
{
433+
cv::eigenNonSymmetric(src, eigenvalues0, noArray());
434+
cv::eigenNonSymmetric(src, eigenvalues, eigenvectors);
435+
}
436+
#if 0
437+
std::cout << "src = " << src << std::endl;
438+
std::cout << "eigenvalues.t() = " << eigenvalues.t() << std::endl;
439+
std::cout << "eigenvectors = " << eigenvectors << std::endl;
440+
#endif
441+
ASSERT_EQ(type, eigenvalues0.type());
442+
ASSERT_EQ(type, eigenvalues.type());
443+
ASSERT_EQ(type, eigenvectors.type());
444+
445+
ASSERT_EQ(src.rows, eigenvalues.rows);
446+
ASSERT_EQ(eigenvalues.rows, eigenvectors.rows);
447+
ASSERT_EQ(src.rows, eigenvectors.cols);
448+
449+
EXPECT_LT(cvtest::norm(eigenvalues, eigenvalues0, NORM_INF), eps);
450+
451+
// check definition: src*eigenvectors.row(i).t() = eigenvalues.at<srcType>(i)*eigenvectors.row(i).t()
452+
for (int i = 0; i < src.rows; i++)
453+
{
454+
EXPECT_NEAR(eigenvalues.at<T>(i), expected_eigenvalues(i), eps) << "i=" << i;
455+
Mat lhs = src*eigenvectors.row(i).t();
456+
Mat rhs = eigenvalues.at<T>(i)*eigenvectors.row(i).t();
457+
EXPECT_LT(cvtest::norm(lhs, rhs, NORM_INF), eps)
458+
<< "i=" << i << " eigenvalue=" << eigenvalues.at<T>(i) << std::endl
459+
<< "lhs=" << lhs.t() << std::endl
460+
<< "rhs=" << rhs.t();
461+
}
462+
}
463+
464+
template<typename T>
465+
static void testEigenSymmetric3x3()
466+
{
467+
/*const*/ T values_[] = {
468+
2, -1, 0,
469+
-1, 2, -1,
470+
0, -1, 2
471+
};
472+
Mat_<T> src(3, 3, values_);
473+
474+
/*const*/ T expected_eigenvalues_[] = { 3.414213562373095f, 2, 0.585786437626905f };
475+
Mat_<T> expected_eigenvalues(3, 1, expected_eigenvalues_);
476+
477+
testEigen(src, expected_eigenvalues);
478+
testEigen(src, expected_eigenvalues, true);
479+
}
480+
TEST(Core_EigenSymmetric, float3x3) { testEigenSymmetric3x3<float>(); }
481+
TEST(Core_EigenSymmetric, double3x3) { testEigenSymmetric3x3<double>(); }
482+
483+
template<typename T>
484+
static void testEigenSymmetric5x5()
485+
{
486+
/*const*/ T values_[5*5] = {
487+
5, -1, 0, 2, 1,
488+
-1, 4, -1, 0, 0,
489+
0, -1, 3, 1, -1,
490+
2, 0, 1, 4, 0,
491+
1, 0, -1, 0, 1
492+
};
493+
Mat_<T> src(5, 5, values_);
494+
495+
/*const*/ T expected_eigenvalues_[] = { 7.028919644935684f, 4.406130784616501f, 3.73626552682258f, 1.438067799899037f, 0.390616243726198f };
496+
Mat_<T> expected_eigenvalues(5, 1, expected_eigenvalues_);
497+
498+
testEigen(src, expected_eigenvalues);
499+
testEigen(src, expected_eigenvalues, true);
500+
}
501+
TEST(Core_EigenSymmetric, float5x5) { testEigenSymmetric5x5<float>(); }
502+
TEST(Core_EigenSymmetric, double5x5) { testEigenSymmetric5x5<double>(); }
503+
504+
505+
template<typename T>
506+
static void testEigen2x2()
507+
{
508+
/*const*/ T values_[] = { 4, 1, 6, 3 };
509+
Mat_<T> src(2, 2, values_);
510+
511+
/*const*/ T expected_eigenvalues_[] = { 6, 1 };
512+
Mat_<T> expected_eigenvalues(2, 1, expected_eigenvalues_);
513+
514+
testEigen(src, expected_eigenvalues);
515+
}
516+
TEST(Core_EigenNonSymmetric, float2x2) { testEigen2x2<float>(); }
517+
TEST(Core_EigenNonSymmetric, double2x2) { testEigen2x2<double>(); }
518+
519+
template<typename T>
520+
static void testEigen3x3()
521+
{
522+
/*const*/ T values_[3*3] = {
523+
3,1,0,
524+
0,3,1,
525+
0,0,3
526+
};
527+
Mat_<T> src(3, 3, values_);
528+
529+
/*const*/ T expected_eigenvalues_[] = { 3, 3, 3 };
530+
Mat_<T> expected_eigenvalues(3, 1, expected_eigenvalues_);
531+
532+
testEigen(src, expected_eigenvalues);
533+
}
534+
TEST(Core_EigenNonSymmetric, float3x3) { testEigen3x3<float>(); }
535+
TEST(Core_EigenNonSymmetric, double3x3) { testEigen3x3<double>(); }

0 commit comments

Comments
 (0)