opencv/tests/cv/src/astereomatching.cpp

   1 /*M///////////////////////////////////////////////////////////////////////////////////////
   2 //
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
   4 //
   5 //  By downloading, copying, installing or using the software you agree to this license.
   6 //  If you do not agree to this license, do not download, install,
   7 //  copy or use the software.
   8 //
   9 //
  10 //                        Intel License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
  14 // Third party copyrights are property of their respective owners.
  15 //
  16 // Redistribution and use in source and binary forms, with or without modification,
  17 // are permitted provided that the following conditions are met:
  18 //
  19 //   * Redistribution's of source code must retain the above copyright notice,
  20 //     this list of conditions and the following disclaimer.
  21 //
  22 //   * Redistribution's in binary form must reproduce the above copyright notice,
  23 //     this list of conditions and the following disclaimer in the documentation
  24 //     and/or other materials provided with the distribution.
  25 //
  26 //   * The name of Intel Corporation may not be used to endorse or promote products
  27 //     derived from this software without specific prior written permission.
  28 //
  29 // This software is provided by the copyright holders and contributors "as is" and
  30 // any express or implied warranties, including, but not limited to, the implied
  31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  32 // In no event shall the Intel Corporation or contributors be liable for any direct,
  33 // indirect, incidental, special, exemplary, or consequential damages
  34 // (including, but not limited to, procurement of substitute goods or services;
  35 // loss of use, data, or profits; or business interruption) however caused
  36 // and on any theory of liability, whether in contract, strict liability,
  37 // or tort (including negligence or otherwise) arising in any way out of
  38 // the use of this software, even if advised of the possibility of such damage.
  39 //
  40 //M*/
  41
  42 /*
  43   This is a regression test for stereo matching algorithms. This test gets some quality metrics
  44   discribed in "A Taxonomy and Evaluation of Dense Two-Frame Stereo Correspondence Algorithms".
  45   Daniel Scharstein, Richard Szeliski
  46 */
  47
  48 #include "cvtest.h"
  49
  50 using namespace std;
  51 using namespace cv;
  52
  53 const float EVAL_BAD_THRESH = 1.f;
  54 const int EVAL_TEXTURELESS_WIDTH = 9;
  55 const float EVAL_TEXTURELESS_THRESH = 2.f;
  56 const float EVAL_DISP_THRESH = 1.f;
  57 const float EVAL_DISP_GAP = 2.f;
  58 const int EVAL_DISCONT_WIDTH = 9;
  59 const int EVAL_IGNORE_BORDER = 10;
  60
  61 const int ERROR_KINDS_COUNT = 6;
  62
  63 //============================== quality measuring functions =================================================
  64
  65 /*
  66   Calculate textureless regions of image: regions where the squared horizontal intensity gradient averaged over
  67   a square window of size=evalTexturelessWidth is below a threshold=evalTexturelessThresh.
  68 */
  69 void computeTexturelessRegions( const Mat& img,  Mat& texturelessMask,
  70              int texturelessWidth = EVAL_TEXTURELESS_WIDTH, float texturelessThresh = EVAL_TEXTURELESS_THRESH )
  71 {
  72     if( img.empty() )
  73         CV_Error( CV_StsBadArg, "img is empty" );
  74     if( img.type() != CV_8UC3 )
  75         CV_Error( CV_StsBadArg, "img.type() must be CV_8UC1" );
  76     Mat dxI; Sobel( img, dxI, CV_32F, 1, 0, 3 );
  77     Mat dxI2; pow( dxI / 8.f, 2, dxI2 );
  78     Mat tmp; cvtColor(dxI2, tmp, CV_BGR2GRAY); dxI2 = tmp;
  79     Mat avgDxI2; boxFilter( dxI2, avgDxI2, CV_32FC1, Size(texturelessWidth,texturelessWidth) );
  80     texturelessMask = avgDxI2 < texturelessThresh;
  81 }
  82
  83 void checkDispMaps( const Mat& disp1, const Mat& disp2 )
  84 {
  85     if( disp1.empty() || disp2.empty() )
  86         CV_Error( CV_StsBadArg, "disp1 or disp2 is empty" );
  87     if( disp1.depth() != CV_8U || disp2.depth() != CV_8U )
  88         CV_Error( CV_StsBadArg, "disp1.depth() and disp2.depth() must be CV_8U" );
  89     Size sz = disp1.size();
  90     if( disp1.cols != sz.width || disp2.rows != sz.height )
  91         CV_Error( CV_StsBadArg, "disp1 and disp2 must have the same size" );
  92 }
  93
  94 /*
  95   Calculate occluded regions of reference image (left image): regions that are occluded in the matching image (right image),
  96   i.e., where the forward-mapped disparity lands at a location with a larger (nearer) disparity.
  97 */
  98 void computeOccludedRegions( const Mat& leftDisp, const Mat& rightDisp, Mat& occludedMask,
  99                              float dispThresh = EVAL_DISP_THRESH )
 100 {
 101     checkDispMaps(leftDisp, rightDisp);
 102     Mat grayLeftDisp; cvtColor(leftDisp, grayLeftDisp, CV_BGR2GRAY);
 103     Mat grayRightDisp; cvtColor(rightDisp, grayRightDisp, CV_BGR2GRAY);
 104
 105     occludedMask.create(leftDisp.size(), CV_8UC1); occludedMask.setTo(Scalar::all(255));
 106     for( int leftY = 0; leftY < leftDisp.rows; leftY++ )
 107         for( int leftX = 0; leftX < leftDisp.cols; leftX++ )
 108         {
 109             int leftDispVal = grayLeftDisp.at<uchar>(leftY, leftX);
 110             int rightX = leftX - leftDispVal, rightY = leftY;
 111             if( rightX < 0 )
 112                 occludedMask.at<uchar>(leftY, leftX) = 0;
 113             else
 114             {
 115                 int rightDispVal = grayRightDisp.at<uchar>(rightY, rightX);
 116                 if( abs(rightDispVal - leftDispVal) > dispThresh )
 117                     occludedMask.at<uchar>(leftY, leftX) = 0;
 118             }
 119         }
 120 }
 121
 122 /*
 123   Calculate depth discontinuty regions: pixels whose neiboring disparities differ by more than
 124   dispGap, dilated by window of width discontWidth.
 125 */
 126 void computeDepthDiscontRegions( const Mat& disp, Mat& depthDiscontMask,
 127                                  float dispGap = EVAL_DISP_GAP, int discontWidth = EVAL_DISCONT_WIDTH )
 128 {
 129     if( disp.empty() )
 130         CV_Error( CV_StsBadArg, "disp is empty" );
 131     if( disp.depth() != CV_8U )
 132         CV_Error( CV_StsBadArg, "disp.depth() must be CV_8U" );
 133     Mat grayDisp; cvtColor(disp, grayDisp, CV_BGR2GRAY);
 134     Mat maxNeighbDisp; dilate(grayDisp, maxNeighbDisp, Mat(3, 3, CV_8UC1, Scalar(1)));
 135     Mat minNeighbDisp; erode(grayDisp, minNeighbDisp, Mat(3, 3, CV_8UC1, Scalar(1)));
 136     depthDiscontMask = max((Mat)(maxNeighbDisp-grayDisp), (Mat)(grayDisp-minNeighbDisp)) > dispGap;
 137     dilate(depthDiscontMask, depthDiscontMask, Mat(discontWidth, discontWidth, CV_8UC1, Scalar(1)));
 138 }
 139
 140 /*
 141   Following functions are for getting evaluation masks excluding a border.
 142 */
 143 Mat borderedAllMask( Size maskSize, int border = EVAL_IGNORE_BORDER )
 144 {
 145     CV_Assert( border >= 0 );
 146     Mat mask(maskSize, CV_8UC1, Scalar(0));
 147     int w = maskSize.width - 2*border, h = maskSize.height - 2*border;
 148     if( w < 0 ||  h < 0 )
 149         mask.setTo(Scalar(0));
 150     else
 151         mask( Rect(Point(border,border),Size(w,h)) ).setTo(Scalar(255));
 152     return mask;
 153 }
 154
 155 void checkMask( const Mat& mask )
 156 {
 157     if( mask.empty() )
 158         CV_Error( CV_StsBadArg, "mask is empty" );
 159     if( mask.type() != CV_8UC1 )
 160         CV_Error( CV_StsBadArg, "mask must have CV_8UC1 type" );
 161 }
 162
 163 void checkMasks( const Mat& mask1, const Mat& mask2 )
 164 {
 165     checkMask(mask1);
 166     checkMask(mask2);
 167     if( mask1.cols != mask2.cols || mask1.rows != mask2.rows )
 168         CV_Error( CV_StsBadArg, "mask1 and mask2 must have the same size" );
 169 }
 170
 171 Mat borderedNoOcclMask( const Mat& occludedMask, int border = EVAL_IGNORE_BORDER )
 172 {
 173     checkMask(occludedMask);
 174     return ~occludedMask & borderedAllMask(occludedMask.size(), border);
 175 }
 176
 177 Mat borderedOcclMask( const Mat& occludedMask, int border = EVAL_IGNORE_BORDER )
 178 {
 179     checkMask(occludedMask);
 180     return occludedMask & borderedAllMask(occludedMask.size(), border);
 181 }
 182
 183 Mat borderedTexturedMask( const Mat& texturelessMask, const Mat& occludedMask, int border = EVAL_IGNORE_BORDER )
 184 {
 185     checkMasks(texturelessMask, occludedMask);
 186     return ~texturelessMask & borderedAllMask(occludedMask.size(), border);
 187 }
 188
 189 Mat borderedTexturelessMask( const Mat& texturelessMask, const Mat& occludedMask, int border = EVAL_IGNORE_BORDER )
 190 {
 191     checkMasks(texturelessMask, occludedMask);
 192     return texturelessMask & borderedAllMask(occludedMask.size(), border);
 193 }
 194
 195 Mat borderedDepthDiscontMask( const Mat& depthDiscontMask, const Mat& occludedMask, int border = EVAL_IGNORE_BORDER )
 196 {
 197     checkMasks(depthDiscontMask, occludedMask);
 198     return depthDiscontMask & borderedAllMask(occludedMask.size(), border);
 199 }
 200
 201 /*
 202   Calculate root-mean-squared error between the computed disparity map (computedDisp) and ground truth map (groundTruthDisp).
 203 */
 204 float dispRMS( const Mat& computedDisp, const Mat& groundTruthDisp, const Mat& mask )
 205 {
 206     checkDispMaps(computedDisp, groundTruthDisp);
 207     Mat grayComputedDisp; cvtColor(computedDisp, grayComputedDisp, CV_BGR2GRAY);
 208     Mat grayGroundTruthDisp; cvtColor(groundTruthDisp, grayGroundTruthDisp, CV_BGR2GRAY);
 209
 210     int pointsCount = mask.empty() ? computedDisp.cols*computedDisp.rows : countNonZero(mask);
 211     return 1.f/sqrt((float)pointsCount) * norm(grayComputedDisp, grayGroundTruthDisp, NORM_L2, mask);
 212 }
 213
 214 /*
 215   Calculate percentage of bad matching pixels.
 216 */
 217 float badMatchPxlsPercentage( const Mat& computedDisp, const Mat& groundTruthDisp, const Mat& mask,
 218                               int badThresh = EVAL_BAD_THRESH )
 219 {
 220     checkDispMaps(computedDisp, groundTruthDisp);
 221     Mat grayComputedDisp; cvtColor(computedDisp, grayComputedDisp, CV_BGR2GRAY);
 222     Mat grayGroundTruthDisp; cvtColor(groundTruthDisp, grayGroundTruthDisp, CV_BGR2GRAY);
 223
 224     Mat badPxlsMap;
 225     absdiff(computedDisp, groundTruthDisp, badPxlsMap);
 226     badPxlsMap = badPxlsMap > badThresh;
 227     int pointsCount = computedDisp.cols*computedDisp.rows;
 228     if( !mask.empty() )
 229     {
 230         badPxlsMap = badPxlsMap & mask;
 231         pointsCount = countNonZero(mask);
 232     }
 233     return 1.f/pointsCount * countNonZero(badPxlsMap);
 234 }
 235
 236 /*
 237   Calculate root-mean-squared errors for six kinds regions:
 238   bordered region, bordered non occluded region, bordered occluded region, bordered textured region,
 239   bordered textureless region, bordered depth discontinuty region.
 240 */
 241 void calcRMSs( const Mat& computedDisp, const Mat& groundTruthDisp,
 242                const Mat& texturelessMask, const Mat& occludedMask, const Mat& depthDiscontMask,
 243                vector<float>& errors )
 244 {
 245     errors.resize(ERROR_KINDS_COUNT);
 246     errors[0] = dispRMS( computedDisp, groundTruthDisp, borderedAllMask(computedDisp.size()) );
 247     errors[1] = dispRMS( computedDisp, groundTruthDisp, borderedNoOcclMask(occludedMask) );
 248     errors[2] = dispRMS( computedDisp, groundTruthDisp, borderedOcclMask(occludedMask) ),
 249     errors[3] = dispRMS( computedDisp, groundTruthDisp, borderedTexturedMask(texturelessMask, occludedMask) ),
 250     errors[4] = dispRMS( computedDisp, groundTruthDisp, borderedTexturelessMask(texturelessMask, occludedMask) ),
 251     errors[5] = dispRMS( computedDisp, groundTruthDisp, borderedDepthDiscontMask(depthDiscontMask, occludedMask) );
 252 }
 253
 254 /*
 255   Calculate percentages of bad matching pixels for six kinds regions:
 256   bordered region, bordered non occluded region, bordered occluded region, bordered textured region,
 257   bordered textureless region, bordered depth discontinuty region.
 258 */
 259 void calcBadMatchPxlsPercentages( const Mat& computedDisp, const Mat& groundTruthDisp,
 260                 const Mat& texturelessMask, const Mat& occludedMask, const Mat& depthDiscontMask,
 261                 vector<float>& errors, int badThresh = EVAL_BAD_THRESH )
 262 {
 263     errors.resize(ERROR_KINDS_COUNT);
 264     errors[0] = badMatchPxlsPercentage( computedDisp, groundTruthDisp, borderedAllMask(computedDisp.size()), badThresh ),
 265     errors[1] = badMatchPxlsPercentage( computedDisp, groundTruthDisp, borderedNoOcclMask(occludedMask), badThresh ),
 266     errors[2] = badMatchPxlsPercentage( computedDisp, groundTruthDisp, borderedOcclMask(occludedMask), badThresh ),
 267     errors[3] = badMatchPxlsPercentage( computedDisp, groundTruthDisp, borderedTexturedMask(texturelessMask, occludedMask), badThresh ),
 268     errors[4] = badMatchPxlsPercentage( computedDisp, groundTruthDisp, borderedTexturelessMask(texturelessMask, occludedMask), badThresh ),
 269     errors[5] = badMatchPxlsPercentage( computedDisp, groundTruthDisp, borderedDepthDiscontMask(depthDiscontMask, occludedMask), badThresh );
 270 }
 271
 272 //===================== regression test for stereo matching algorithms ==============================
 273
 274 const string RESULT_DIR = "test_results";
 275 const string LEFT_IMG_NAME = "im2.ppm";
 276 const string RIGHT_IMG_NAME = "im6.ppm";
 277 const string TRUE_LEFT_DISP_NAME = "disp2.pgm";
 278 const string TRUE_RIGHT_DISP_NAME = "disp6.pgm";
 279
 280 string DATASETS_NAMES[] = { "barn2", "bull", "cones", "poster", "sawtooth", "teddy", "venus" };
 281 string ERROR_PREFIXES[] = { "borderedAll",
 282                             "borderedNoOccl",
 283                             "borderedOccl",
 284                             "borderedTextured",
 285                             "borderedTextureless",
 286                             "borderedDepthDiscont" }; // size of ERROR_KINDS_COUNT
 287
 288 const string RMS_STR = "RMS";
 289 const string BAD_PXLS_PERCENTAGE_STR = "badPxlsPercentage";
 290
 291 class CV_StereoMatchingTest : public CvTest
 292 {
 293 public:
 294     CV_StereoMatchingTest(const char* testName, const char* testFuncs) :
 295             CvTest( testName, testFuncs ) {}
 296 protected:
 297     // assumed that left image is a reference image
 298     virtual void runStereoMatchingFunc( const Mat& leftImg, const Mat& rigthImg,
 299                                         Mat& leftDisp, Mat& rightDisp ) = 0;
 300
 301     void writeErrors( FileStorage fs, const string& errName, const vector<float>& errors );
 302     void readErrors( FileNode fn, const string& errName, vector<float>& errors );
 303     int compareErrors( const vector<float>& calcErrors, const vector<float>& validErrors,
 304                        const vector<float>& eps, const string& errName, const string& datasetName );
 305     int processStereoMatchingResults( FileStorage& fs, const string& datasetName, bool isWrite,
 306                   const Mat& leftImg, const Mat& rightImg,
 307                   const Mat& trueLeftDisp, const Mat& trueRightDisp,
 308                   const Mat& leftDisp, const Mat& rightDisp );
 309     void run( int );
 310
 311     string resultFilename;
 312 };
 313
 314 void CV_StereoMatchingTest::writeErrors( FileStorage fs, const string& errName, const vector<float>& errors )
 315 {
 316     assert( (int)errors.size() == ERROR_KINDS_COUNT );
 317     vector<float>::const_iterator it = errors.begin();
 318     for( int i = 0; i < ERROR_KINDS_COUNT; i++, ++it )
 319         fs << ERROR_PREFIXES[i] + errName << *it;
 320 }
 321
 322 void CV_StereoMatchingTest::readErrors( FileNode fn, const string& errName, vector<float>& errors )
 323 {
 324     errors.resize( ERROR_KINDS_COUNT );
 325     vector<float>::iterator it = errors.begin();
 326     for( int i = 0; i < ERROR_KINDS_COUNT; i++, ++it )
 327         fn[ERROR_PREFIXES[i]+errName] >> *it;
 328 }
 329
 330 int CV_StereoMatchingTest::compareErrors( const vector<float>& calcErrors, const vector<float>& validErrors,
 331                    const vector<float>& eps, const string& errName, const string& datasetName )
 332 {
 333     assert( (int)calcErrors.size() == ERROR_KINDS_COUNT );
 334     assert( (int)validErrors.size() == ERROR_KINDS_COUNT );
 335     assert( (int)eps.size() == ERROR_KINDS_COUNT );
 336     vector<float>::const_iterator calcIt = calcErrors.begin(),
 337                                   validIt = validErrors.begin(),
 338                                   epsIt = eps.begin();
 339     for( int i = 0; i < ERROR_KINDS_COUNT; i++, ++calcIt, ++validIt, ++epsIt )
 340         if( fabs(*calcIt - *validIt) > *epsIt )
 341         {
 342             ts->printf( CvTS::LOG, "bad accuracy of %s on dataset %s",
 343                         string(ERROR_PREFIXES[i]+errName).c_str(), datasetName.c_str() );
 344             return CvTS::FAIL_BAD_ACCURACY;
 345         }
 346     return CvTS::OK;
 347 }
 348
 349 int CV_StereoMatchingTest::processStereoMatchingResults( FileStorage& fs, const string& datasetName, bool isWrite,
 350               const Mat& leftImg, const Mat& rightImg,
 351               const Mat& trueLeftDisp, const Mat& trueRightDisp,
 352               const Mat& leftDisp, const Mat& rightDisp )
 353 {
 354     // rightDisp is not used in current test virsion
 355     int code = CvTS::OK;
 356     assert( fs.isOpened() );
 357     assert( !datasetName.empty() );
 358
 359     // commpute three kinds masks
 360     Mat texturelessMask; computeTexturelessRegions( leftImg, texturelessMask );
 361     Mat occludedMask; computeOccludedRegions( trueLeftDisp, trueRightDisp, occludedMask );
 362     Mat depthDiscontMask; computeDepthDiscontRegions( trueLeftDisp, depthDiscontMask );
 363     assert( !texturelessMask.empty() && !occludedMask.empty() && !depthDiscontMask.empty() );
 364
 365     // compute errors
 366     vector<float> rmss, badPxlsPercentages;
 367     calcRMSs( leftDisp, trueLeftDisp, texturelessMask, occludedMask, depthDiscontMask, rmss );
 368     calcBadMatchPxlsPercentages( leftDisp, trueLeftDisp,
 369               texturelessMask, occludedMask, depthDiscontMask, badPxlsPercentages );
 370
 371     if( isWrite )
 372     {
 373         fs << datasetName << "{";
 374         writeErrors( fs, RMS_STR, rmss );
 375         writeErrors( fs, BAD_PXLS_PERCENTAGE_STR, badPxlsPercentages );
 376         fs << "}"; // datasetName
 377     }
 378     else // compare
 379     {
 380         FileNode fn = fs.getFirstTopLevelNode()[datasetName];
 381         vector<float> validRmss, validBadPxlsPercentages;
 382
 383         readErrors( fn, RMS_STR, validRmss );
 384         int tempCode = compareErrors( rmss, validRmss,
 385                       vector<float>(ERROR_KINDS_COUNT, 0.01f), RMS_STR, datasetName );
 386         code = tempCode==CvTS::OK ? code : tempCode;
 387
 388         readErrors( fn, BAD_PXLS_PERCENTAGE_STR, validBadPxlsPercentages );
 389         tempCode = compareErrors( badPxlsPercentages, validBadPxlsPercentages,
 390                       vector<float>(ERROR_KINDS_COUNT, 0.01f), BAD_PXLS_PERCENTAGE_STR, datasetName );
 391         code = tempCode==CvTS::OK ? code : tempCode;
 392     }
 393     return code;
 394 }
 395
 396 void CV_StereoMatchingTest::run(int)
 397 {
 398     int code = CvTS::OK;
 399     if( resultFilename.empty() )
 400     {
 401         ts->printf( CvTS::LOG, "resultFilename is empty" );
 402         ts->set_failed_test_info( CvTS::FAIL_BAD_ARG_CHECK );
 403         return;
 404     }
 405
 406     string fullResultFilename = string(ts->get_data_path()) + "stereomatching/" + RESULT_DIR + "/" + resultFilename;
 407     bool isWrite = true; // write or compare results
 408     FileStorage fs( fullResultFilename, FileStorage::READ );
 409     if( fs.isOpened() )
 410         isWrite = false;
 411     else
 412     {
 413         fs.open( fullResultFilename, FileStorage::WRITE );
 414         if( !fs.isOpened() )
 415         {
 416             ts->printf( CvTS::LOG, "file named %s can not be read or written", fullResultFilename.c_str() );
 417             code = CvTS::FAIL_BAD_ARG_CHECK;
 418             return;
 419         }
 420         fs << "stereo_matching" << "{";
 421     }
 422
 423     int datasetsCount = sizeof(DATASETS_NAMES)/sizeof(DATASETS_NAMES[0]);
 424     for( int dsi = 0; dsi < datasetsCount; dsi++)
 425     {
 426         string dataPath = string(ts->get_data_path()) + "stereomatching/datasets/" + DATASETS_NAMES[dsi] + "/";
 427         Mat leftImg = imread(dataPath + LEFT_IMG_NAME);
 428         Mat rightImg = imread(dataPath + RIGHT_IMG_NAME);
 429         Mat trueLeftDisp = imread(dataPath + TRUE_LEFT_DISP_NAME);
 430         Mat trueRightDisp = imread(dataPath + TRUE_RIGHT_DISP_NAME);
 431         if( leftImg.empty() || rightImg.empty() || trueLeftDisp.empty() || trueRightDisp.empty() )
 432         {
 433             ts->printf( CvTS::LOG, "images or true disparities of dataset %s can not be read", DATASETS_NAMES[dsi].c_str() );
 434             code = CvTS::FAIL_INVALID_TEST_DATA;
 435             continue;
 436         }
 437
 438         Mat leftDisp, rightDisp;
 439         runStereoMatchingFunc( leftImg, rightImg, leftDisp, rightDisp );
 440
 441         int tempCode = processStereoMatchingResults( fs, DATASETS_NAMES[dsi], isWrite,
 442                    leftImg, rightImg, trueLeftDisp, trueRightDisp, leftDisp, rightDisp);
 443         code = tempCode==CvTS::OK ? code : tempCode;
 444     }
 445
 446     if( isWrite )
 447         fs << "}"; // "stereo_matching"
 448
 449     ts->set_failed_test_info( code );
 450 }
 451
 452 //CV_StereoMatchingTest stereoMatchingTest;