Skip to content

Commit ce20efb

Browse files
alalekalalek
committed
Merge pull request opencv#9804 from woodychow:optimize_cveigen
2 parents b0bce60 + c0b6061 commit ce20efb

File tree

4 files changed

+292
-292
lines changed

4 files changed

+292
-292
lines changed

‎modules/core/misc/java/test/CoreTest.java

Lines changed: 51 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -394,26 +394,36 @@ public void testDivideMatMatMatDoubleInt() {
394394
}
395395

396396
public void testEigen() {
397-
Mat src = new Mat(3, 3, CvType.CV_32FC1, new Scalar(2.0));
397+
Mat src = new Mat(3, 3, CvType.CV_32FC1) {
398+
{
399+
put(0, 0, 2, 0, 0);
400+
put(1, 0, 0, 6, 0);
401+
put(2, 0, 0, 0, 4);
402+
}
403+
};
398404
Mat eigenVals = new Mat();
399405
Mat eigenVecs = new Mat();
400406

401407
Core.eigen(src, eigenVals, eigenVecs);
402408

403409
Mat expectedEigenVals = new Mat(3, 1, CvType.CV_32FC1) {
404410
{
405-
put(0, 0, 6, 0, 0);
406-
}
407-
};
408-
Mat expectedEigenVecs = new Mat(3, 3, CvType.CV_32FC1) {
409-
{
410-
put(0, 0, 0.57735026, 0.57735026, 0.57735032);
411-
put(1, 0, 0.70710677, -0.70710677, 0);
412-
put(2, 0, -0.40824831, -0.40824831, 0.81649661);
411+
put(0, 0, 6, 4, 2);
413412
}
414413
};
415414
assertMatEqual(eigenVals, expectedEigenVals, EPS);
416-
assertMatEqual(eigenVecs, expectedEigenVecs, EPS);
415+
416+
// check by definition
417+
double eps = 1e-3;
418+
for(int i = 0; i < 3; i++)
419+
{
420+
Mat vec = eigenVecs.row(i).t();
421+
Mat lhs = new Mat(3, 1, CvType.CV_32FC1);
422+
Core.gemm(src, vec, 1.0, new Mat(), 1.0, lhs);
423+
Mat rhs = new Mat(3, 1, CvType.CV_32FC1);
424+
Core.gemm(vec, eigenVals.row(i), 1.0, new Mat(), 1.0, rhs);
425+
assertMatEqual(lhs, rhs, eps);
426+
}
417427
}
418428

419429
public void testExp() {
@@ -1326,7 +1336,8 @@ public void testPCAComputeMatMatMat() {
13261336
Mat vectors = new Mat();
13271337

13281338
Core.PCACompute(data, mean, vectors);
1329-
1339+
//System.out.println(mean.dump());
1340+
//System.out.println(vectors.dump());
13301341
Mat mean_truth = new Mat(1, 4, CvType.CV_32F) {
13311342
{
13321343
put(0, 0, 2, 4, 4, 8);
@@ -1338,7 +1349,21 @@ public void testPCAComputeMatMatMat() {
13381349
}
13391350
};
13401351
assertMatEqual(mean_truth, mean, EPS);
1341-
assertMatEqual(vectors_truth, vectors, EPS);
1352+
1353+
// eigenvectors are normalized (length = 1),
1354+
// but direction is unknown (v and -v are both eigen vectors)
1355+
// so this direct check doesn't work:
1356+
// assertMatEqual(vectors_truth, vectors, EPS);
1357+
for(int i = 0; i < 3; i++)
1358+
{
1359+
Mat vec0 = vectors_truth.row(i);
1360+
Mat vec1 = vectors.row(i);
1361+
Mat vec1_ = new Mat();
1362+
Core.subtract(new Mat(1, 4, CvType.CV_32F, new Scalar(0)), vec1, vec1_);
1363+
double scale1 = Core.norm(vec0, vec1);
1364+
double scale2 = Core.norm(vec0, vec1_);
1365+
assertTrue(Math.min(scale1, scale2) < EPS);
1366+
}
13421367
}
13431368

13441369
public void testPCAComputeMatMatMatInt() {
@@ -1365,7 +1390,20 @@ public void testPCAComputeMatMatMatInt() {
13651390
}
13661391
};
13671392
assertMatEqual(mean_truth, mean, EPS);
1368-
assertMatEqual(vectors_truth, vectors, EPS);
1393+
// eigenvectors are normalized (length = 1),
1394+
// but direction is unknown (v and -v are both eigen vectors)
1395+
// so this direct check doesn't work:
1396+
// assertMatEqual(vectors_truth, vectors, EPS);
1397+
for(int i = 0; i < 1; i++)
1398+
{
1399+
Mat vec0 = vectors_truth.row(i);
1400+
Mat vec1 = vectors.row(i);
1401+
Mat vec1_ = new Mat();
1402+
Core.subtract(new Mat(1, 4, CvType.CV_32F, new Scalar(0)), vec1, vec1_);
1403+
double scale1 = Core.norm(vec0, vec1);
1404+
double scale2 = Core.norm(vec0, vec1_);
1405+
assertTrue(Math.min(scale1, scale2) < EPS);
1406+
}
13691407
}
13701408

13711409
public void testPCAProject() {

‎modules/core/src/lapack.cpp

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -43,6 +43,12 @@
4343
#include "precomp.hpp"
4444
#include <limits>
4545

46+
#ifdef HAVE_EIGEN
47+
#include <Eigen/Core>
48+
#include <Eigen/Eigenvalues>
49+
#include "opencv2/core/eigen.hpp"
50+
#endif
51+
4652
#if defined _M_IX86 && defined _MSC_VER && _MSC_VER < 1700
4753
#pragma float_control(precise, on)
4854
#endif
@@ -1396,6 +1402,47 @@ bool cv::eigen( InputArray _src, OutputArray _evals, OutputArray _evects )
13961402
v = _evects.getMat();
13971403
}
13981404

1405+
#ifdef HAVE_EIGEN
1406+
const bool evecNeeded = _evects.needed();
1407+
const int esOptions = evecNeeded ? Eigen::ComputeEigenvectors : Eigen::EigenvaluesOnly;
1408+
_evals.create(n, 1, type);
1409+
cv::Mat evals = _evals.getMat();
1410+
if ( type == CV_64F )
1411+
{
1412+
Eigen::MatrixXd src_eig, zeros_eig;
1413+
cv::cv2eigen(src, src_eig);
1414+
1415+
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es;
1416+
es.compute(src_eig, esOptions);
1417+
if ( es.info() == Eigen::Success )
1418+
{
1419+
cv::eigen2cv(es.eigenvalues().reverse().eval(), evals);
1420+
if ( evecNeeded )
1421+
{
1422+
cv::Mat evects = _evects.getMat();
1423+
cv::eigen2cv(es.eigenvectors().rowwise().reverse().transpose().eval(), v);
1424+
}
1425+
return true;
1426+
}
1427+
} else { // CV_32F
1428+
Eigen::MatrixXf src_eig, zeros_eig;
1429+
cv::cv2eigen(src, src_eig);
1430+
1431+
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXf> es;
1432+
es.compute(src_eig, esOptions);
1433+
if ( es.info() == Eigen::Success )
1434+
{
1435+
cv::eigen2cv(es.eigenvalues().reverse().eval(), evals);
1436+
if ( evecNeeded )
1437+
{
1438+
cv::eigen2cv(es.eigenvectors().rowwise().reverse().transpose().eval(), v);
1439+
}
1440+
return true;
1441+
}
1442+
}
1443+
return false;
1444+
#else
1445+
13991446
size_t elemSize = src.elemSize(), astep = alignSize(n*elemSize, 16);
14001447
AutoBuffer<uchar> buf(n*astep + n*5*elemSize + 32);
14011448
uchar* ptr = alignPtr((uchar*)buf, 16);
@@ -1408,6 +1455,7 @@ bool cv::eigen( InputArray _src, OutputArray _evals, OutputArray _evects )
14081455

14091456
w.copyTo(_evals);
14101457
return ok;
1458+
#endif
14111459
}
14121460

14131461
namespace cv

‎modules/core/test/test_eigen.cpp

Lines changed: 22 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -59,7 +59,7 @@ using namespace std;
5959
#define MESSAGE_ERROR_DIFF_1 "Accuracy of eigen values computing less than required."
6060
#define MESSAGE_ERROR_DIFF_2 "Accuracy of eigen vectors computing less than required."
6161
#define MESSAGE_ERROR_ORTHO "Matrix of eigen vectors is not orthogonal."
62-
#define MESSAGE_ERROR_ORDER "Eigen values are not sorted in ascending order."
62+
#define MESSAGE_ERROR_ORDER "Eigen values are not sorted in descending order."
6363

6464
const int COUNT_NORM_TYPES = 3;
6565
const int NORM_TYPE[COUNT_NORM_TYPES] = {cv::NORM_L1, cv::NORM_L2, cv::NORM_INF};
@@ -164,8 +164,8 @@ void Core_EigenTest_32::run(int) { check_full(CV_32FC1); }
164164
void Core_EigenTest_64::run(int) { check_full(CV_64FC1); }
165165

166166
Core_EigenTest::Core_EigenTest()
167-
: eps_val_32(1e-3f), eps_vec_32(12e-3f),
168-
eps_val_64(1e-4f), eps_vec_64(1e-3f), ntests(100) {}
167+
: eps_val_32(1e-3f), eps_vec_32(1e-3f),
168+
eps_val_64(1e-4f), eps_vec_64(1e-4f), ntests(100) {}
169169
Core_EigenTest::~Core_EigenTest() {}
170170

171171
bool Core_EigenTest::check_pair_count(const cv::Mat& src, const cv::Mat& evalues, int low_index, int high_index)
@@ -234,7 +234,7 @@ bool Core_EigenTest::check_orthogonality(const cv::Mat& U)
234234

235235
for (int i = 0; i < COUNT_NORM_TYPES; ++i)
236236
{
237-
double diff = cvtest::norm(UUt, E, NORM_TYPE[i]);
237+
double diff = cvtest::norm(UUt, E, NORM_TYPE[i] | cv::NORM_RELATIVE);
238238
if (diff > eps_vec)
239239
{
240240
std::cout << endl; std::cout << "Checking orthogonality of matrix " << U << ": ";
@@ -257,7 +257,7 @@ bool Core_EigenTest::check_pairs_order(const cv::Mat& eigen_values)
257257
if (!(eigen_values.at<float>(i, 0) > eigen_values.at<float>(i+1, 0)))
258258
{
259259
std::cout << endl; std::cout << "Checking order of eigen values vector " << eigen_values << "..." << endl;
260-
std::cout << "Pair of indexes with non ascending of eigen values: (" << i << ", " << i+1 << ")." << endl;
260+
std::cout << "Pair of indexes with non descending of eigen values: (" << i << ", " << i+1 << ")." << endl;
261261
std::cout << endl;
262262
CV_Error(CORE_EIGEN_ERROR_ORDER, MESSAGE_ERROR_ORDER);
263263
return false;
@@ -272,9 +272,9 @@ bool Core_EigenTest::check_pairs_order(const cv::Mat& eigen_values)
272272
if (!(eigen_values.at<double>(i, 0) > eigen_values.at<double>(i+1, 0)))
273273
{
274274
std::cout << endl; std::cout << "Checking order of eigen values vector " << eigen_values << "..." << endl;
275-
std::cout << "Pair of indexes with non ascending of eigen values: (" << i << ", " << i+1 << ")." << endl;
275+
std::cout << "Pair of indexes with non descending of eigen values: (" << i << ", " << i+1 << ")." << endl;
276276
std::cout << endl;
277-
CV_Error(CORE_EIGEN_ERROR_ORDER, "Eigen values are not sorted in ascending order.");
277+
CV_Error(CORE_EIGEN_ERROR_ORDER, "Eigen values are not sorted in descending order.");
278278
return false;
279279
}
280280

@@ -307,43 +307,28 @@ bool Core_EigenTest::test_pairs(const cv::Mat& src)
307307

308308
cv::Mat eigen_vectors_t; cv::transpose(eigen_vectors, eigen_vectors_t);
309309

310-
cv::Mat src_evec(src.rows, src.cols, type);
311-
src_evec = src*eigen_vectors_t;
310+
// Check:
311+
// src * eigenvector = eigenval * eigenvector
312+
cv::Mat lhs(src.rows, src.cols, type);
313+
cv::Mat rhs(src.rows, src.cols, type);
312314

313-
cv::Mat eval_evec(src.rows, src.cols, type);
315+
lhs = src*eigen_vectors_t;
314316

315-
switch (type)
317+
for (int i = 0; i < src.cols; ++i)
316318
{
317-
case CV_32FC1:
318-
{
319-
for (int i = 0; i < src.cols; ++i)
320-
{
321-
cv::Mat tmp = eigen_values.at<float>(i, 0) * eigen_vectors_t.col(i);
322-
for (int j = 0; j < src.rows; ++j) eval_evec.at<float>(j, i) = tmp.at<float>(j, 0);
323-
}
324-
325-
break;
326-
}
327-
328-
case CV_64FC1:
319+
double eigenval = 0;
320+
switch (type)
329321
{
330-
for (int i = 0; i < src.cols; ++i)
331-
{
332-
cv::Mat tmp = eigen_values.at<double>(i, 0) * eigen_vectors_t.col(i);
333-
for (int j = 0; j < src.rows; ++j) eval_evec.at<double>(j, i) = tmp.at<double>(j, 0);
334-
}
335-
336-
break;
322+
case CV_32FC1: eigenval = eigen_values.at<float>(i, 0); break;
323+
case CV_64FC1: eigenval = eigen_values.at<double>(i, 0); break;
337324
}
338-
339-
default:;
325+
cv::Mat rhs_v = eigenval * eigen_vectors_t.col(i);
326+
rhs_v.copyTo(rhs.col(i));
340327
}
341328

342-
cv::Mat disparity = src_evec - eval_evec;
343-
344329
for (int i = 0; i < COUNT_NORM_TYPES; ++i)
345330
{
346-
double diff = cvtest::norm(disparity, NORM_TYPE[i]);
331+
double diff = cvtest::norm(lhs, rhs, NORM_TYPE[i] | cv::NORM_RELATIVE);
347332
if (diff > eps_vec)
348333
{
349334
std::cout << endl; std::cout << "Checking accuracy of eigen vectors computing for matrix " << src << ": ";
@@ -372,7 +357,7 @@ bool Core_EigenTest::test_values(const cv::Mat& src)
372357

373358
for (int i = 0; i < COUNT_NORM_TYPES; ++i)
374359
{
375-
double diff = cvtest::norm(eigen_values_1, eigen_values_2, NORM_TYPE[i]);
360+
double diff = cvtest::norm(eigen_values_1, eigen_values_2, NORM_TYPE[i] | cv::NORM_RELATIVE);
376361
if (diff > eps_val)
377362
{
378363
std::cout << endl; std::cout << "Checking accuracy of eigen values computing for matrix " << src << ": ";
@@ -419,7 +404,7 @@ static void testEigen(const Mat_<T>& src, const Mat_<T>& expected_eigenvalues, b
419404
SCOPED_TRACE(runSymmetric ? "cv::eigen" : "cv::eigenNonSymmetric");
420405

421406
int type = traits::Type<T>::value;
422-
const T eps = 1e-6f;
407+
const T eps = src.type() == CV_32F ? 1e-4f : 1e-6f;
423408

424409
Mat eigenvalues, eigenvectors, eigenvalues0;
425410

0 commit comments

Comments
 (0)