diff --git a/modules/core/misc/java/test/CoreTest.java b/modules/core/misc/java/test/CoreTest.java index 3187e7cd2d..bd05407e20 100644 --- a/modules/core/misc/java/test/CoreTest.java +++ b/modules/core/misc/java/test/CoreTest.java @@ -394,7 +394,13 @@ public class CoreTest extends OpenCVTestCase { } public void testEigen() { - Mat src = new Mat(3, 3, CvType.CV_32FC1, new Scalar(2.0)); + Mat src = new Mat(3, 3, CvType.CV_32FC1) { + { + put(0, 0, 2, 0, 0); + put(1, 0, 0, 6, 0); + put(2, 0, 0, 0, 4); + } + }; Mat eigenVals = new Mat(); Mat eigenVecs = new Mat(); @@ -402,18 +408,22 @@ public class CoreTest extends OpenCVTestCase { Mat expectedEigenVals = new Mat(3, 1, CvType.CV_32FC1) { { - put(0, 0, 6, 0, 0); - } - }; - Mat expectedEigenVecs = new Mat(3, 3, CvType.CV_32FC1) { - { - put(0, 0, 0.57735026, 0.57735026, 0.57735032); - put(1, 0, 0.70710677, -0.70710677, 0); - put(2, 0, -0.40824831, -0.40824831, 0.81649661); + put(0, 0, 6, 4, 2); } }; assertMatEqual(eigenVals, expectedEigenVals, EPS); - assertMatEqual(eigenVecs, expectedEigenVecs, EPS); + + // check by definition + double eps = 1e-3; + for(int i = 0; i < 3; i++) + { + Mat vec = eigenVecs.row(i).t(); + Mat lhs = new Mat(3, 1, CvType.CV_32FC1); + Core.gemm(src, vec, 1.0, new Mat(), 1.0, lhs); + Mat rhs = new Mat(3, 1, CvType.CV_32FC1); + Core.gemm(vec, eigenVals.row(i), 1.0, new Mat(), 1.0, rhs); + assertMatEqual(lhs, rhs, eps); + } } public void testExp() { @@ -1326,7 +1336,8 @@ public class CoreTest extends OpenCVTestCase { Mat vectors = new Mat(); Core.PCACompute(data, mean, vectors); - + //System.out.println(mean.dump()); + //System.out.println(vectors.dump()); Mat mean_truth = new Mat(1, 4, CvType.CV_32F) { { put(0, 0, 2, 4, 4, 8); @@ -1338,7 +1349,21 @@ public class CoreTest extends OpenCVTestCase { } }; assertMatEqual(mean_truth, mean, EPS); - assertMatEqual(vectors_truth, vectors, EPS); + + // eigenvectors are normalized (length = 1), + // but direction is unknown (v and -v are both eigen vectors) + // so this direct check doesn't work: + // assertMatEqual(vectors_truth, vectors, EPS); + for(int i = 0; i < 3; i++) + { + Mat vec0 = vectors_truth.row(i); + Mat vec1 = vectors.row(i); + Mat vec1_ = new Mat(); + Core.subtract(new Mat(1, 4, CvType.CV_32F, new Scalar(0)), vec1, vec1_); + double scale1 = Core.norm(vec0, vec1); + double scale2 = Core.norm(vec0, vec1_); + assertTrue(Math.min(scale1, scale2) < EPS); + } } public void testPCAComputeMatMatMatInt() { @@ -1365,7 +1390,20 @@ public class CoreTest extends OpenCVTestCase { } }; assertMatEqual(mean_truth, mean, EPS); - assertMatEqual(vectors_truth, vectors, EPS); + // eigenvectors are normalized (length = 1), + // but direction is unknown (v and -v are both eigen vectors) + // so this direct check doesn't work: + // assertMatEqual(vectors_truth, vectors, EPS); + for(int i = 0; i < 1; i++) + { + Mat vec0 = vectors_truth.row(i); + Mat vec1 = vectors.row(i); + Mat vec1_ = new Mat(); + Core.subtract(new Mat(1, 4, CvType.CV_32F, new Scalar(0)), vec1, vec1_); + double scale1 = Core.norm(vec0, vec1); + double scale2 = Core.norm(vec0, vec1_); + assertTrue(Math.min(scale1, scale2) < EPS); + } } public void testPCAProject() { diff --git a/modules/core/src/lapack.cpp b/modules/core/src/lapack.cpp index 800eb0ec87..8fdbddad6e 100644 --- a/modules/core/src/lapack.cpp +++ b/modules/core/src/lapack.cpp @@ -43,6 +43,12 @@ #include "precomp.hpp" #include +#ifdef HAVE_EIGEN +#include +#include +#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 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 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 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 index 44f4a729c8..b57fe94c18 100644 --- 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}; @@ -164,8 +164,8 @@ void Core_EigenTest_32::run(int) { check_full(CV_32FC1); } void Core_EigenTest_64::run(int) { check_full(CV_64FC1); } Core_EigenTest::Core_EigenTest() -: eps_val_32(1e-3f), eps_vec_32(12e-3f), - eps_val_64(1e-4f), eps_vec_64(1e-3f), ntests(100) {} +: eps_val_32(1e-3f), eps_vec_32(1e-3f), + eps_val_64(1e-4f), eps_vec_64(1e-4f), ntests(100) {} Core_EigenTest::~Core_EigenTest() {} 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) for (int i = 0; i < COUNT_NORM_TYPES; ++i) { - double diff = cvtest::norm(UUt, E, NORM_TYPE[i]); + double diff = cvtest::norm(UUt, E, NORM_TYPE[i] | cv::NORM_RELATIVE); if (diff > eps_vec) { 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) if (!(eigen_values.at(i, 0) > eigen_values.at(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(i, 0) > eigen_values.at(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; } @@ -307,43 +307,28 @@ bool Core_EigenTest::test_pairs(const cv::Mat& src) cv::Mat eigen_vectors_t; cv::transpose(eigen_vectors, eigen_vectors_t); - cv::Mat src_evec(src.rows, src.cols, type); - src_evec = src*eigen_vectors_t; + // Check: + // src * eigenvector = eigenval * eigenvector + cv::Mat lhs(src.rows, src.cols, type); + cv::Mat rhs(src.rows, src.cols, type); - cv::Mat eval_evec(src.rows, src.cols, type); + lhs = src*eigen_vectors_t; - switch (type) + for (int i = 0; i < src.cols; ++i) { - case CV_32FC1: - { - for (int i = 0; i < src.cols; ++i) - { - cv::Mat tmp = eigen_values.at(i, 0) * eigen_vectors_t.col(i); - for (int j = 0; j < src.rows; ++j) eval_evec.at(j, i) = tmp.at(j, 0); - } - - break; - } - - case CV_64FC1: + double eigenval = 0; + switch (type) { - for (int i = 0; i < src.cols; ++i) - { - cv::Mat tmp = eigen_values.at(i, 0) * eigen_vectors_t.col(i); - for (int j = 0; j < src.rows; ++j) eval_evec.at(j, i) = tmp.at(j, 0); - } - - break; + case CV_32FC1: eigenval = eigen_values.at(i, 0); break; + case CV_64FC1: eigenval = eigen_values.at(i, 0); break; } - - default:; + cv::Mat rhs_v = eigenval * eigen_vectors_t.col(i); + rhs_v.copyTo(rhs.col(i)); } - cv::Mat disparity = src_evec - eval_evec; - for (int i = 0; i < COUNT_NORM_TYPES; ++i) { - double diff = cvtest::norm(disparity, NORM_TYPE[i]); + double diff = cvtest::norm(lhs, rhs, NORM_TYPE[i] | cv::NORM_RELATIVE); if (diff > eps_vec) { 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) for (int i = 0; i < COUNT_NORM_TYPES; ++i) { - double diff = cvtest::norm(eigen_values_1, eigen_values_2, NORM_TYPE[i]); + double diff = cvtest::norm(eigen_values_1, eigen_values_2, NORM_TYPE[i] | cv::NORM_RELATIVE); if (diff > eps_val) { std::cout << endl; std::cout << "Checking accuracy of eigen values computing for matrix " << src << ": "; @@ -419,7 +404,7 @@ static void testEigen(const Mat_& src, const Mat_& expected_eigenvalues, b SCOPED_TRACE(runSymmetric ? "cv::eigen" : "cv::eigenNonSymmetric"); int type = traits::Type::value; - const T eps = 1e-6f; + const T eps = src.type() == CV_32F ? 1e-4f : 1e-6f; Mat eigenvalues, eigenvectors, eigenvalues0; diff --git a/modules/core/test/test_mat.cpp b/modules/core/test/test_mat.cpp index 8d04831a23..b4abc52323 100644 --- a/modules/core/test/test_mat.cpp +++ b/modules/core/test/test_mat.cpp @@ -286,258 +286,188 @@ void Core_ReduceTest::run( int ) #define CHECK_C -class Core_PCATest : public cvtest::BaseTest +TEST(Core_PCA, accuracy) { -public: - Core_PCATest() {} -protected: - void run(int) - { - const Size sz(200, 500); - - double diffPrjEps, diffBackPrjEps, - prjEps, backPrjEps, - evalEps, evecEps; - int maxComponents = 100; - double retainedVariance = 0.95; - Mat rPoints(sz, CV_32FC1), rTestPoints(sz, CV_32FC1); - RNG& rng = ts->get_rng(); - - rng.fill( rPoints, RNG::UNIFORM, Scalar::all(0.0), Scalar::all(1.0) ); - rng.fill( rTestPoints, RNG::UNIFORM, Scalar::all(0.0), Scalar::all(1.0) ); - - PCA rPCA( rPoints, Mat(), CV_PCA_DATA_AS_ROW, maxComponents ), cPCA; - - // 1. check C++ PCA & ROW - Mat rPrjTestPoints = rPCA.project( rTestPoints ); - Mat rBackPrjTestPoints = rPCA.backProject( rPrjTestPoints ); - - Mat avg(1, sz.width, CV_32FC1 ); - cv::reduce( rPoints, avg, 0, CV_REDUCE_AVG ); - Mat Q = rPoints - repeat( avg, rPoints.rows, 1 ), Qt = Q.t(), eval, evec; - Q = Qt * Q; - Q = Q /(float)rPoints.rows; - - eigen( Q, eval, evec ); - /*SVD svd(Q); - evec = svd.vt; - eval = svd.w;*/ - - Mat subEval( maxComponents, 1, eval.type(), eval.ptr() ), - subEvec( maxComponents, evec.cols, evec.type(), evec.ptr() ); - - #ifdef CHECK_C - Mat prjTestPoints, backPrjTestPoints, cPoints = rPoints.t(), cTestPoints = rTestPoints.t(); - CvMat _points, _testPoints, _avg, _eval, _evec, _prjTestPoints, _backPrjTestPoints; - #endif - - // check eigen() - double eigenEps = 1e-6; - double err; - for(int i = 0; i < Q.rows; i++ ) - { - Mat v = evec.row(i).t(); - Mat Qv = Q * v; - - Mat lv = eval.at(i,0) * v; - err = cvtest::norm( Qv, lv, NORM_L2 ); - if( err > eigenEps ) - { - ts->printf( cvtest::TS::LOG, "bad accuracy of eigen(); err = %f\n", err ); - ts->set_failed_test_info( cvtest::TS::FAIL_BAD_ACCURACY ); - return; - } - } - // check pca eigenvalues - evalEps = 1e-6, evecEps = 1e-3; - err = cvtest::norm( rPCA.eigenvalues, subEval, NORM_L2 ); - if( err > evalEps ) - { - ts->printf( cvtest::TS::LOG, "pca.eigenvalues is incorrect (CV_PCA_DATA_AS_ROW); err = %f\n", err ); - ts->set_failed_test_info( cvtest::TS::FAIL_BAD_ACCURACY ); - return; - } - // check pca eigenvectors - for(int i = 0; i < subEvec.rows; i++) - { - Mat r0 = rPCA.eigenvectors.row(i); - Mat r1 = subEvec.row(i); - err = cvtest::norm( r0, r1, CV_L2 ); - if( err > evecEps ) - { - r1 *= -1; - double err2 = cvtest::norm(r0, r1, CV_L2); - if( err2 > evecEps ) - { - Mat tmp; - absdiff(rPCA.eigenvectors, subEvec, tmp); - double mval = 0; Point mloc; - minMaxLoc(tmp, 0, &mval, 0, &mloc); - - ts->printf( cvtest::TS::LOG, "pca.eigenvectors is incorrect (CV_PCA_DATA_AS_ROW); err = %f\n", err ); - ts->printf( cvtest::TS::LOG, "max diff is %g at (i=%d, j=%d) (%g vs %g)\n", - mval, mloc.y, mloc.x, rPCA.eigenvectors.at(mloc.y, mloc.x), - subEvec.at(mloc.y, mloc.x)); - ts->set_failed_test_info( cvtest::TS::FAIL_BAD_ACCURACY ); - return; - } - } - } - - prjEps = 1.265, backPrjEps = 1.265; - for( int i = 0; i < rTestPoints.rows; i++ ) - { - // check pca project - Mat subEvec_t = subEvec.t(); - Mat prj = rTestPoints.row(i) - avg; prj *= subEvec_t; - err = cvtest::norm(rPrjTestPoints.row(i), prj, CV_RELATIVE_L2); - if( err > prjEps ) - { - ts->printf( cvtest::TS::LOG, "bad accuracy of project() (CV_PCA_DATA_AS_ROW); err = %f\n", err ); - ts->set_failed_test_info( cvtest::TS::FAIL_BAD_ACCURACY ); - return; - } - // check pca backProject - Mat backPrj = rPrjTestPoints.row(i) * subEvec + avg; - err = cvtest::norm( rBackPrjTestPoints.row(i), backPrj, CV_RELATIVE_L2 ); - if( err > backPrjEps ) - { - ts->printf( cvtest::TS::LOG, "bad accuracy of backProject() (CV_PCA_DATA_AS_ROW); err = %f\n", err ); - ts->set_failed_test_info( cvtest::TS::FAIL_BAD_ACCURACY ); - return; - } - } - - // 2. check C++ PCA & COL - cPCA( rPoints.t(), Mat(), CV_PCA_DATA_AS_COL, maxComponents ); - diffPrjEps = 1, diffBackPrjEps = 1; - Mat ocvPrjTestPoints = cPCA.project(rTestPoints.t()); - err = cvtest::norm(cv::abs(ocvPrjTestPoints), cv::abs(rPrjTestPoints.t()), CV_RELATIVE_L2 ); - if( err > diffPrjEps ) - { - ts->printf( cvtest::TS::LOG, "bad accuracy of project() (CV_PCA_DATA_AS_COL); err = %f\n", err ); - ts->set_failed_test_info( cvtest::TS::FAIL_BAD_ACCURACY ); - return; - } - err = cvtest::norm(cPCA.backProject(ocvPrjTestPoints), rBackPrjTestPoints.t(), CV_RELATIVE_L2 ); - if( err > diffBackPrjEps ) - { - ts->printf( cvtest::TS::LOG, "bad accuracy of backProject() (CV_PCA_DATA_AS_COL); err = %f\n", err ); - ts->set_failed_test_info( cvtest::TS::FAIL_BAD_ACCURACY ); - return; - } - - // 3. check C++ PCA w/retainedVariance - cPCA( rPoints.t(), Mat(), CV_PCA_DATA_AS_COL, retainedVariance ); - diffPrjEps = 1, diffBackPrjEps = 1; - Mat rvPrjTestPoints = cPCA.project(rTestPoints.t()); - - if( cPCA.eigenvectors.rows > maxComponents) - err = cvtest::norm(cv::abs(rvPrjTestPoints.rowRange(0,maxComponents)), cv::abs(rPrjTestPoints.t()), CV_RELATIVE_L2 ); - else - err = cvtest::norm(cv::abs(rvPrjTestPoints), cv::abs(rPrjTestPoints.colRange(0,cPCA.eigenvectors.rows).t()), CV_RELATIVE_L2 ); + const Size sz(200, 500); + + double diffPrjEps, diffBackPrjEps, + prjEps, backPrjEps, + evalEps, evecEps; + int maxComponents = 100; + double retainedVariance = 0.95; + Mat rPoints(sz, CV_32FC1), rTestPoints(sz, CV_32FC1); + RNG rng(12345); + + rng.fill( rPoints, RNG::UNIFORM, Scalar::all(0.0), Scalar::all(1.0) ); + rng.fill( rTestPoints, RNG::UNIFORM, Scalar::all(0.0), Scalar::all(1.0) ); + + PCA rPCA( rPoints, Mat(), CV_PCA_DATA_AS_ROW, maxComponents ), cPCA; + + // 1. check C++ PCA & ROW + Mat rPrjTestPoints = rPCA.project( rTestPoints ); + Mat rBackPrjTestPoints = rPCA.backProject( rPrjTestPoints ); + + Mat avg(1, sz.width, CV_32FC1 ); + cv::reduce( rPoints, avg, 0, CV_REDUCE_AVG ); + Mat Q = rPoints - repeat( avg, rPoints.rows, 1 ), Qt = Q.t(), eval, evec; + Q = Qt * Q; + Q = Q /(float)rPoints.rows; + + eigen( Q, eval, evec ); + /*SVD svd(Q); + evec = svd.vt; + eval = svd.w;*/ + + Mat subEval( maxComponents, 1, eval.type(), eval.ptr() ), + subEvec( maxComponents, evec.cols, evec.type(), evec.ptr() ); + +#ifdef CHECK_C + Mat prjTestPoints, backPrjTestPoints, cPoints = rPoints.t(), cTestPoints = rTestPoints.t(); + CvMat _points, _testPoints, _avg, _eval, _evec, _prjTestPoints, _backPrjTestPoints; +#endif - if( err > diffPrjEps ) - { - ts->printf( cvtest::TS::LOG, "bad accuracy of project() (CV_PCA_DATA_AS_COL); retainedVariance=0.95; err = %f\n", err ); - ts->set_failed_test_info( cvtest::TS::FAIL_BAD_ACCURACY ); - return; - } - err = cvtest::norm(cPCA.backProject(rvPrjTestPoints), rBackPrjTestPoints.t(), CV_RELATIVE_L2 ); - if( err > diffBackPrjEps ) - { - ts->printf( cvtest::TS::LOG, "bad accuracy of backProject() (CV_PCA_DATA_AS_COL); retainedVariance=0.95; err = %f\n", err ); - ts->set_failed_test_info( cvtest::TS::FAIL_BAD_ACCURACY ); - return; - } + // check eigen() + double eigenEps = 1e-4; + double err; + for(int i = 0; i < Q.rows; i++ ) + { + Mat v = evec.row(i).t(); + Mat Qv = Q * v; - #ifdef CHECK_C - // 4. check C PCA & ROW - _points = rPoints; - _testPoints = rTestPoints; - _avg = avg; - _eval = eval; - _evec = evec; - prjTestPoints.create(rTestPoints.rows, maxComponents, rTestPoints.type() ); - backPrjTestPoints.create(rPoints.size(), rPoints.type() ); - _prjTestPoints = prjTestPoints; - _backPrjTestPoints = backPrjTestPoints; - - cvCalcPCA( &_points, &_avg, &_eval, &_evec, CV_PCA_DATA_AS_ROW ); - cvProjectPCA( &_testPoints, &_avg, &_evec, &_prjTestPoints ); - cvBackProjectPCA( &_prjTestPoints, &_avg, &_evec, &_backPrjTestPoints ); - - err = cvtest::norm(prjTestPoints, rPrjTestPoints, CV_RELATIVE_L2); - if( err > diffPrjEps ) - { - ts->printf( cvtest::TS::LOG, "bad accuracy of cvProjectPCA() (CV_PCA_DATA_AS_ROW); err = %f\n", err ); - ts->set_failed_test_info( cvtest::TS::FAIL_BAD_ACCURACY ); - return; - } - err = cvtest::norm(backPrjTestPoints, rBackPrjTestPoints, CV_RELATIVE_L2); - if( err > diffBackPrjEps ) + Mat lv = eval.at(i,0) * v; + err = cvtest::norm(Qv, lv, NORM_L2 | NORM_RELATIVE); + EXPECT_LE(err, eigenEps) << "bad accuracy of eigen(); i = " << i; + } + // check pca eigenvalues + evalEps = 1e-5, evecEps = 5e-3; + err = cvtest::norm(rPCA.eigenvalues, subEval, NORM_L2 | NORM_RELATIVE); + EXPECT_LE(err , evalEps) << "pca.eigenvalues is incorrect (CV_PCA_DATA_AS_ROW)"; + // check pca eigenvectors + for(int i = 0; i < subEvec.rows; i++) + { + Mat r0 = rPCA.eigenvectors.row(i); + Mat r1 = subEvec.row(i); + // eigenvectors have normalized length, but both directions v and -v are valid + double err1 = cvtest::norm(r0, r1, NORM_L2 | NORM_RELATIVE); + double err2 = cvtest::norm(r0, -r1, NORM_L2 | NORM_RELATIVE); + err = std::min(err1, err2); + if (err > evecEps) { - ts->printf( cvtest::TS::LOG, "bad accuracy of cvBackProjectPCA() (CV_PCA_DATA_AS_ROW); err = %f\n", err ); - ts->set_failed_test_info( cvtest::TS::FAIL_BAD_ACCURACY ); - return; + Mat tmp; + absdiff(rPCA.eigenvectors, subEvec, tmp); + double mval = 0; Point mloc; + minMaxLoc(tmp, 0, &mval, 0, &mloc); + + EXPECT_LE(err, evecEps) << "pca.eigenvectors is incorrect (CV_PCA_DATA_AS_ROW) at " << i << " " + << cv::format("max diff is %g at (i=%d, j=%d) (%g vs %g)\n", + mval, mloc.y, mloc.x, rPCA.eigenvectors.at(mloc.y, mloc.x), + subEvec.at(mloc.y, mloc.x)) + << "r0=" << r0 << std::endl + << "r1=" << r1 << std::endl + << "err1=" << err1 << " err2=" << err2 + ; } + } - // 5. check C PCA & COL - _points = cPoints; - _testPoints = cTestPoints; - avg = avg.t(); _avg = avg; - eval = eval.t(); _eval = eval; - evec = evec.t(); _evec = evec; - prjTestPoints = prjTestPoints.t(); _prjTestPoints = prjTestPoints; - backPrjTestPoints = backPrjTestPoints.t(); _backPrjTestPoints = backPrjTestPoints; - - cvCalcPCA( &_points, &_avg, &_eval, &_evec, CV_PCA_DATA_AS_COL ); - cvProjectPCA( &_testPoints, &_avg, &_evec, &_prjTestPoints ); - cvBackProjectPCA( &_prjTestPoints, &_avg, &_evec, &_backPrjTestPoints ); - - err = cvtest::norm(cv::abs(prjTestPoints), cv::abs(rPrjTestPoints.t()), CV_RELATIVE_L2 ); - if( err > diffPrjEps ) - { - ts->printf( cvtest::TS::LOG, "bad accuracy of cvProjectPCA() (CV_PCA_DATA_AS_COL); err = %f\n", err ); - ts->set_failed_test_info( cvtest::TS::FAIL_BAD_ACCURACY ); - return; - } - err = cvtest::norm(backPrjTestPoints, rBackPrjTestPoints.t(), CV_RELATIVE_L2); - if( err > diffBackPrjEps ) - { - ts->printf( cvtest::TS::LOG, "bad accuracy of cvBackProjectPCA() (CV_PCA_DATA_AS_COL); err = %f\n", err ); - ts->set_failed_test_info( cvtest::TS::FAIL_BAD_ACCURACY ); - return; - } - #endif - // Test read and write - FileStorage fs( "PCA_store.yml", FileStorage::WRITE ); - rPCA.write( fs ); - fs.release(); - - PCA lPCA; - fs.open( "PCA_store.yml", FileStorage::READ ); - lPCA.read( fs.root() ); - err = cvtest::norm( rPCA.eigenvectors, lPCA.eigenvectors, CV_RELATIVE_L2 ); - if( err > 0 ) - { - ts->printf( cvtest::TS::LOG, "bad accuracy of write/load functions (YML); err = %f\n", err ); - ts->set_failed_test_info( cvtest::TS::FAIL_BAD_ACCURACY ); - } - err = cvtest::norm( rPCA.eigenvalues, lPCA.eigenvalues, CV_RELATIVE_L2 ); - if( err > 0 ) + prjEps = 1.265, backPrjEps = 1.265; + for( int i = 0; i < rTestPoints.rows; i++ ) + { + // check pca project + Mat subEvec_t = subEvec.t(); + Mat prj = rTestPoints.row(i) - avg; prj *= subEvec_t; + err = cvtest::norm(rPrjTestPoints.row(i), prj, NORM_L2 | NORM_RELATIVE); + if (err < prjEps) { - ts->printf( cvtest::TS::LOG, "bad accuracy of write/load functions (YML); err = %f\n", err ); - ts->set_failed_test_info( cvtest::TS::FAIL_BAD_ACCURACY ); + EXPECT_LE(err, prjEps) << "bad accuracy of project() (CV_PCA_DATA_AS_ROW)"; + continue; } - err = cvtest::norm( rPCA.mean, lPCA.mean, CV_RELATIVE_L2 ); - if( err > 0 ) + // check pca backProject + Mat backPrj = rPrjTestPoints.row(i) * subEvec + avg; + err = cvtest::norm(rBackPrjTestPoints.row(i), backPrj, NORM_L2 | NORM_RELATIVE); + if (err > backPrjEps) { - ts->printf( cvtest::TS::LOG, "bad accuracy of write/load functions (YML); err = %f\n", err ); - ts->set_failed_test_info( cvtest::TS::FAIL_BAD_ACCURACY ); + EXPECT_LE(err, backPrjEps) << "bad accuracy of backProject() (CV_PCA_DATA_AS_ROW)"; + continue; } } -}; + + // 2. check C++ PCA & COL + cPCA( rPoints.t(), Mat(), CV_PCA_DATA_AS_COL, maxComponents ); + diffPrjEps = 1, diffBackPrjEps = 1; + Mat ocvPrjTestPoints = cPCA.project(rTestPoints.t()); + err = cvtest::norm(cv::abs(ocvPrjTestPoints), cv::abs(rPrjTestPoints.t()), NORM_L2 | NORM_RELATIVE); + ASSERT_LE(err, diffPrjEps) << "bad accuracy of project() (CV_PCA_DATA_AS_COL)"; + err = cvtest::norm(cPCA.backProject(ocvPrjTestPoints), rBackPrjTestPoints.t(), NORM_L2 | NORM_RELATIVE); + ASSERT_LE(err, diffBackPrjEps) << "bad accuracy of backProject() (CV_PCA_DATA_AS_COL)"; + + // 3. check C++ PCA w/retainedVariance + cPCA( rPoints.t(), Mat(), CV_PCA_DATA_AS_COL, retainedVariance ); + diffPrjEps = 1, diffBackPrjEps = 1; + Mat rvPrjTestPoints = cPCA.project(rTestPoints.t()); + + if( cPCA.eigenvectors.rows > maxComponents) + err = cvtest::norm(cv::abs(rvPrjTestPoints.rowRange(0,maxComponents)), cv::abs(rPrjTestPoints.t()), NORM_L2 | NORM_RELATIVE); + else + err = cvtest::norm(cv::abs(rvPrjTestPoints), cv::abs(rPrjTestPoints.colRange(0,cPCA.eigenvectors.rows).t()), NORM_L2 | NORM_RELATIVE); + + ASSERT_LE(err, diffPrjEps) << "bad accuracy of project() (CV_PCA_DATA_AS_COL); retainedVariance=" << retainedVariance; + err = cvtest::norm(cPCA.backProject(rvPrjTestPoints), rBackPrjTestPoints.t(), NORM_L2 | NORM_RELATIVE); + ASSERT_LE(err, diffBackPrjEps) << "bad accuracy of backProject() (CV_PCA_DATA_AS_COL); retainedVariance=" << retainedVariance; + +#ifdef CHECK_C + // 4. check C PCA & ROW + _points = rPoints; + _testPoints = rTestPoints; + _avg = avg; + _eval = eval; + _evec = evec; + prjTestPoints.create(rTestPoints.rows, maxComponents, rTestPoints.type() ); + backPrjTestPoints.create(rPoints.size(), rPoints.type() ); + _prjTestPoints = prjTestPoints; + _backPrjTestPoints = backPrjTestPoints; + + cvCalcPCA( &_points, &_avg, &_eval, &_evec, CV_PCA_DATA_AS_ROW ); + cvProjectPCA( &_testPoints, &_avg, &_evec, &_prjTestPoints ); + cvBackProjectPCA( &_prjTestPoints, &_avg, &_evec, &_backPrjTestPoints ); + + err = cvtest::norm(prjTestPoints, rPrjTestPoints, NORM_L2 | NORM_RELATIVE); + ASSERT_LE(err, diffPrjEps) << "bad accuracy of cvProjectPCA() (CV_PCA_DATA_AS_ROW)"; + err = cvtest::norm(backPrjTestPoints, rBackPrjTestPoints, NORM_L2 | NORM_RELATIVE); + ASSERT_LE(err, diffBackPrjEps) << "bad accuracy of cvBackProjectPCA() (CV_PCA_DATA_AS_ROW)"; + + // 5. check C PCA & COL + _points = cPoints; + _testPoints = cTestPoints; + avg = avg.t(); _avg = avg; + eval = eval.t(); _eval = eval; + evec = evec.t(); _evec = evec; + prjTestPoints = prjTestPoints.t(); _prjTestPoints = prjTestPoints; + backPrjTestPoints = backPrjTestPoints.t(); _backPrjTestPoints = backPrjTestPoints; + + cvCalcPCA( &_points, &_avg, &_eval, &_evec, CV_PCA_DATA_AS_COL ); + cvProjectPCA( &_testPoints, &_avg, &_evec, &_prjTestPoints ); + cvBackProjectPCA( &_prjTestPoints, &_avg, &_evec, &_backPrjTestPoints ); + + err = cvtest::norm(cv::abs(prjTestPoints), cv::abs(rPrjTestPoints.t()), NORM_L2 | NORM_RELATIVE); + ASSERT_LE(err, diffPrjEps) << "bad accuracy of cvProjectPCA() (CV_PCA_DATA_AS_COL)"; + err = cvtest::norm(backPrjTestPoints, rBackPrjTestPoints.t(), NORM_L2 | NORM_RELATIVE); + ASSERT_LE(err, diffBackPrjEps) << "bad accuracy of cvBackProjectPCA() (CV_PCA_DATA_AS_COL)"; +#endif + // Test read and write + FileStorage fs( "PCA_store.yml", FileStorage::WRITE ); + rPCA.write( fs ); + fs.release(); + + PCA lPCA; + fs.open( "PCA_store.yml", FileStorage::READ ); + lPCA.read( fs.root() ); + err = cvtest::norm(rPCA.eigenvectors, lPCA.eigenvectors, NORM_L2 | NORM_RELATIVE); + EXPECT_LE(err, 0) << "bad accuracy of write/load functions (YML)"; + err = cvtest::norm(rPCA.eigenvalues, lPCA.eigenvalues, NORM_L2 | NORM_RELATIVE); + EXPECT_LE(err, 0) << "bad accuracy of write/load functions (YML)"; + err = cvtest::norm(rPCA.mean, lPCA.mean, NORM_L2 | NORM_RELATIVE); + EXPECT_LE(err, 0) << "bad accuracy of write/load functions (YML)"; +} class Core_ArrayOpTest : public cvtest::BaseTest { @@ -1227,7 +1157,6 @@ protected: } }; -TEST(Core_PCA, accuracy) { Core_PCATest test; test.safe_run(); } TEST(Core_Reduce, accuracy) { Core_ReduceTest test; test.safe_run(); } TEST(Core_Array, basic_operations) { Core_ArrayOpTest test; test.safe_run(); }