opencv/tests/cxcore/src/apca.cpp

   1 /*M///////////////////////////////////////////////////////////////////////////////////////\r
   2 //\r
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.\r
   4 //\r
   5 //  By downloading, copying, installing or using the software you agree to this license.\r
   6 //  If you do not agree to this license, do not download, install,\r
   7 //  copy or use the software.\r
   8 //\r
   9 //\r
  10 //                        Intel License Agreement\r
  11 //                For Open Source Computer Vision Library\r
  12 //\r
  13 // Copyright (C) 2000, Intel Corporation, all rights reserved.\r
  14 // Third party copyrights are property of their respective owners.\r
  15 //\r
  16 // Redistribution and use in source and binary forms, with or without modification,\r
  17 // are permitted provided that the following conditions are met:\r
  18 //\r
  19 //   * Redistribution's of source code must retain the above copyright notice,\r
  20 //     this list of conditions and the following disclaimer.\r
  21 //\r
  22 //   * Redistribution's in binary form must reproduce the above copyright notice,\r
  23 //     this list of conditions and the following disclaimer in the documentation\r
  24 //     and/or other materials provided with the distribution.\r
  25 //\r
  26 //   * The name of Intel Corporation may not be used to endorse or promote products\r
  27 //     derived from this software without specific prior written permission.\r
  28 //\r
  29 // This software is provided by the copyright holders and contributors "as is" and\r
  30 // any express or implied warranties, including, but not limited to, the implied\r
  31 // warranties of merchantability and fitness for a particular purpose are disclaimed.\r
  32 // In no event shall the Intel Corporation or contributors be liable for any direct,\r
  33 // indirect, incidental, special, exemplary, or consequential damages\r
  34 // (including, but not limited to, procurement of substitute goods or services;\r
  35 // loss of use, data, or profits; or business interruption) however caused\r
  36 // and on any theory of liability, whether in contract, strict liability,\r
  37 // or tort (including negligence or otherwise) arising in any way out of\r
  38 // the use of this software, even if advised of the possibility of such damage.\r
  39 //\r
  40 //M*/\r
  41 \r
  42 #include "cxcoretest.h"\r
  43 \r
  44 using namespace cv;\r
  45 \r
  46 #define CHECK_C\r
  47 \r
  48 Size sz( 200, 500);\r
  49 \r
  50 class CV_PCATest : public CvTest\r
  51 {\r
  52 public:\r
  53     CV_PCATest() : CvTest( "pca", "PCA funcs" ) {}\r
  54 protected:\r
  55     void run( int);\r
  56 };\r
  57 \r
  58 #if 0\r
  59 \r
  60 void CV_PCATest::run( int )\r
  61 {\r
  62     int code = CvTS::OK, err;\r
  63     int maxComponents = 1;\r
  64     Mat points( 1000, 3, CV_32FC1);\r
  65 \r
  66         RNG rng = *ts->get_rng(); // get ts->rng seed\r
  67         rng.fill( points, RNG::NORMAL, Scalar::all(0.0), Scalar::all(1.0) );\r
  68 \r
  69     float mp[] = { 3.0f, 3.0f, 3.0f }, cp[] = { 0.5f, 0.0f, 0.0f,\r
  70                                                 0.0f, 1.0f, 0.0f,\r
  71                                                 0.0f, 0.0f, 0.3f };\r
  72     Mat mean( 1, 3, CV_32FC1, mp ),\r
  73         cov( 3, 3, CV_32FC1, cp );\r
  74     for( int i = 0; i < points.rows; i++ )\r
  75     {\r
  76         Mat r(1, points.cols, CV_32FC1, points.ptr<float>(i));\r
  77         r =  r * cov + mean; \r
  78     }\r
  79 \r
  80     PCA pca( points, Mat(), CV_PCA_DATA_AS_ROW, maxComponents );\r
  81 \r
  82     // check project\r
  83     Mat prjPoints = pca.project( points );\r
  84     err = 0;\r
  85     for( int i = 0; i < prjPoints.rows; i++ )\r
  86     {\r
  87         float val = prjPoints.at<float>(i,0);\r
  88         if( val > 3.0f || val < -3.0f )\r
  89             err++;\r
  90     }\r
  91         float projectErr = 0.02f;\r
  92         if( (float)err > prjPoints.rows * projectErr )\r
  93     {\r
  94         ts->printf( CvTS::LOG, "bad accuracy of project() (real = %f, permissible = %f)",\r
  95                         (float)err/(float)prjPoints.rows, projectErr );\r
  96         code = CvTS::FAIL_BAD_ACCURACY;\r
  97     }\r
  98 \r
  99     // check backProject\r
 100     Mat points1 = pca.backProject( prjPoints );\r
 101     err = 0;\r
 102         for( int i = 0; i < points.rows; i++ ) \r
 103         {\r
 104                 if( fabs(points1.at<float>(i,0) - mean.at<float>(0,0)) > 0.15 ||\r
 105             fabs(points1.at<float>(i,1) - points.at<float>(i,1)) > 0.05 ||\r
 106             fabs(points1.at<float>(i,2) - mean.at<float>(0,2)) > 0.15 )\r
 107             err++;\r
 108         }\r
 109         float backProjectErr = 0.05f;\r
 110         if( (float)err > prjPoints.rows*backProjectErr )\r
 111     {\r
 112         ts->printf( CvTS::LOG, "bad accuracy of backProject() (real = %f, permissible = %f)",\r
 113                         (float)err/(float)prjPoints.rows, backProjectErr );\r
 114         code = CvTS::FAIL_BAD_ACCURACY;\r
 115     }\r
 116 \r
 117         CvRNG *oldRng = ts->get_rng(); // set ts->rng seed\r
 118         *oldRng = rng.state;\r
 119 \r
 120     ts->set_failed_test_info( code );\r
 121 }\r
 122 #else\r
 123 void CV_PCATest::run( int )\r
 124 {\r
 125         int code = CvTS::OK;\r
 126         \r
 127         double diffPrjEps, diffBackPrjEps,\r
 128                    prjEps, backPrjEps,\r
 129                    evalEps, evecEps;\r
 130         int maxComponents = 100;\r
 131         Mat rPoints(sz, CV_32FC1), rTestPoints(sz, CV_32FC1);\r
 132         RNG rng = *ts->get_rng(); \r
 133 \r
 134         rng.fill( rPoints, RNG::UNIFORM, Scalar::all(0.0), Scalar::all(1.0) );\r
 135         rng.fill( rTestPoints, RNG::UNIFORM, Scalar::all(0.0), Scalar::all(1.0) );\r
 136 \r
 137         PCA rPCA( rPoints, Mat(), CV_PCA_DATA_AS_ROW, maxComponents ), cPCA;\r
 138 \r
 139         // 1. check C++ PCA & ROW\r
 140         Mat rPrjTestPoints = rPCA.project( rTestPoints );\r
 141         Mat rBackPrjTestPoints = rPCA.backProject( rPrjTestPoints );\r
 142 \r
 143         Mat avg(1, sz.width, CV_32FC1 );\r
 144         reduce( rPoints, avg, 0, CV_REDUCE_AVG );\r
 145         Mat Q = rPoints - repeat( avg, rPoints.rows, 1 ), Qt = Q.t(), eval, evec;\r
 146         Q = Qt * Q;\r
 147         Q = Q /(float)rPoints.rows;\r
 148 \r
 149         eigen( Q, eval, evec );\r
 150         /*SVD svd(Q);\r
 151         evec = svd.vt;\r
 152         eval = svd.w;*/\r
 153 \r
 154         Mat subEval( maxComponents, 1, eval.type(), eval.data ),\r
 155                 subEvec( maxComponents, evec.cols, evec.type(), evec.data );\r
 156 \r
 157 #ifdef CHECK_C\r
 158         Mat prjTestPoints, backPrjTestPoints, cPoints = rPoints.t(), cTestPoints = rTestPoints.t();\r
 159         CvMat _points, _testPoints, _avg, _eval, _evec, _prjTestPoints, _backPrjTestPoints;\r
 160 #endif\r
 161 \r
 162         // check eigen()\r
 163         double eigenEps = 1e-6;\r
 164         double err;\r
 165         for(int i = 0; i < Q.rows; i++ )\r
 166         {\r
 167                 Mat v = evec.row(i).t();\r
 168                 Mat Qv = Q * v;\r
 169 \r
 170                 Mat lv = eval.at<float>(i,0) * v;\r
 171                 err = norm( Qv, lv );\r
 172                 if( err > eigenEps )\r
 173                 {\r
 174                         ts->printf( CvTS::LOG, "bad accuracy of eigen(); err = %f\n", err );\r
 175                         code = CvTS::FAIL_BAD_ACCURACY;\r
 176                         goto exit_func;\r
 177                 }\r
 178         }\r
 179         // check pca eigenvalues\r
 180         evalEps = 1e-6, evecEps = 1;\r
 181         err = norm( rPCA.eigenvalues, subEval );\r
 182         if( err > evalEps )\r
 183         {\r
 184                 ts->printf( CvTS::LOG, "pca.eigenvalues is incorrect (CV_PCA_DATA_AS_ROW); err = %f\n", err );\r
 185                 code = CvTS::FAIL_BAD_ACCURACY;\r
 186                 goto exit_func;\r
 187         }\r
 188         // check pca eigenvectors\r
 189         err = norm( rPCA.eigenvectors, subEvec, CV_RELATIVE_L2 );\r
 190         if( err > evecEps )\r
 191         {\r
 192                 ts->printf( CvTS::LOG, "pca.eigenvectors is incorrect (CV_PCA_DATA_AS_ROW); err = %f\n", err );\r
 193                 code = CvTS::FAIL_BAD_ACCURACY;\r
 194                 goto exit_func;\r
 195         }\r
 196         \r
 197     prjEps = 1.265, backPrjEps = 1.265;\r
 198         for( int i = 0; i < rTestPoints.rows; i++ )\r
 199         {\r
 200                 // check pca project\r
 201                 Mat subEvec_t = subEvec.t();\r
 202                 Mat prj = rTestPoints.row(i) - avg; prj *= subEvec_t;\r
 203                 err = norm(rPrjTestPoints.row(i), prj, CV_RELATIVE_L2);\r
 204                 if( err > prjEps )\r
 205                 {\r
 206                         ts->printf( CvTS::LOG, "bad accuracy of project() (CV_PCA_DATA_AS_ROW); err = %f\n", err );\r
 207                         code = CvTS::FAIL_BAD_ACCURACY;\r
 208                         goto exit_func;\r
 209                 }\r
 210                 // check pca backProject\r
 211                 Mat backPrj = rPrjTestPoints.row(i) * subEvec + avg;\r
 212                 err = norm( rBackPrjTestPoints.row(i), backPrj, CV_RELATIVE_L2 );\r
 213                 if( err > backPrjEps )\r
 214                 {\r
 215                         ts->printf( CvTS::LOG, "bad accuracy of backProject() (CV_PCA_DATA_AS_ROW); err = %f\n", err );\r
 216                         code = CvTS::FAIL_BAD_ACCURACY;\r
 217                         goto exit_func;\r
 218                 }\r
 219         }\r
 220 \r
 221         // 2. check C++ PCA & COL\r
 222         cPCA( rPoints.t(), Mat(), CV_PCA_DATA_AS_COL, maxComponents );\r
 223         diffPrjEps = 1, diffBackPrjEps = 1;\r
 224         err = norm(cPCA.project(rTestPoints.t()), rPrjTestPoints.t(), CV_RELATIVE_L2 );\r
 225         if( err > diffPrjEps )\r
 226         {\r
 227                 ts->printf( CvTS::LOG, "bad accuracy of project() (CV_PCA_DATA_AS_COL); err = %f\n", err );\r
 228                 code = CvTS::FAIL_BAD_ACCURACY;\r
 229                 goto exit_func;\r
 230         }\r
 231         err = norm(cPCA.backProject(rPrjTestPoints.t()), rBackPrjTestPoints.t(), CV_RELATIVE_L2 );\r
 232         if( err > diffBackPrjEps )\r
 233         {\r
 234                 ts->printf( CvTS::LOG, "bad accuracy of backProject() (CV_PCA_DATA_AS_COL); err = %f\n", err );\r
 235                 code = CvTS::FAIL_BAD_ACCURACY;\r
 236                 goto exit_func;\r
 237         }\r
 238 \r
 239 #ifdef CHECK_C\r
 240         // 3. check C PCA & ROW\r
 241         _points = rPoints;\r
 242         _testPoints = rTestPoints;\r
 243         _avg = avg;\r
 244         _eval = eval;\r
 245         _evec = evec;\r
 246         prjTestPoints.create(rTestPoints.rows, maxComponents, rTestPoints.type() );\r
 247         backPrjTestPoints.create(rPoints.size(), rPoints.type() );\r
 248         _prjTestPoints = prjTestPoints;\r
 249         _backPrjTestPoints = backPrjTestPoints;\r
 250 \r
 251         cvCalcPCA( &_points, &_avg, &_eval, &_evec, CV_PCA_DATA_AS_ROW );\r
 252         cvProjectPCA( &_testPoints, &_avg, &_evec, &_prjTestPoints );\r
 253         cvBackProjectPCA( &_prjTestPoints, &_avg, &_evec, &_backPrjTestPoints );\r
 254 \r
 255         err = norm(prjTestPoints, rPrjTestPoints, CV_RELATIVE_L2);\r
 256         if( err > diffPrjEps )\r
 257         {\r
 258                 ts->printf( CvTS::LOG, "bad accuracy of cvProjectPCA() (CV_PCA_DATA_AS_ROW); err = %f\n", err );\r
 259                 code = CvTS::FAIL_BAD_ACCURACY;\r
 260                 goto exit_func;\r
 261         }\r
 262         err = norm(backPrjTestPoints, rBackPrjTestPoints, CV_RELATIVE_L2);\r
 263         if( err > diffBackPrjEps )\r
 264         {\r
 265                 ts->printf( CvTS::LOG, "bad accuracy of cvBackProjectPCA() (CV_PCA_DATA_AS_ROW); err = %f\n", err );\r
 266                 code = CvTS::FAIL_BAD_ACCURACY;\r
 267                 goto exit_func;\r
 268         }\r
 269 \r
 270         // 3. check C PCA & COL\r
 271         _points = cPoints;\r
 272         _testPoints = cTestPoints;\r
 273         avg = avg.t(); _avg = avg;\r
 274         eval = eval.t(); _eval = eval;\r
 275         evec = evec.t(); _evec = evec;\r
 276     prjTestPoints = prjTestPoints.t(); _prjTestPoints = prjTestPoints;\r
 277     backPrjTestPoints = backPrjTestPoints.t(); _backPrjTestPoints = backPrjTestPoints;\r
 278 \r
 279         cvCalcPCA( &_points, &_avg, &_eval, &_evec, CV_PCA_DATA_AS_COL );\r
 280         cvProjectPCA( &_testPoints, &_avg, &_evec, &_prjTestPoints );\r
 281         cvBackProjectPCA( &_prjTestPoints, &_avg, &_evec, &_backPrjTestPoints );\r
 282 \r
 283         err = norm(prjTestPoints, rPrjTestPoints.t(), CV_RELATIVE_L2 );\r
 284         if( err > diffPrjEps )\r
 285         {\r
 286                 ts->printf( CvTS::LOG, "bad accuracy of cvProjectPCA() (CV_PCA_DATA_AS_COL); err = %f\n", err );\r
 287                 code = CvTS::FAIL_BAD_ACCURACY;\r
 288                 goto exit_func;\r
 289         }\r
 290         err = norm(backPrjTestPoints, rBackPrjTestPoints.t(), CV_RELATIVE_L2);\r
 291         if( err > diffBackPrjEps )\r
 292         {\r
 293                 ts->printf( CvTS::LOG, "bad accuracy of cvBackProjectPCA() (CV_PCA_DATA_AS_COL); err = %f\n", err );\r
 294                 code = CvTS::FAIL_BAD_ACCURACY;\r
 295                 goto exit_func;\r
 296         }\r
 297 #endif\r
 298 \r
 299 exit_func:\r
 300 \r
 301         CvRNG* _rng = ts->get_rng(); \r
 302         *_rng = rng.state;\r
 303         ts->set_failed_test_info( code );\r
 304 }\r
 305 \r
 306 #endif\r
 307 \r
 308 CV_PCATest pca_test;\r