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 // License Agreement 11 // For Open Source Computer Vision Library 12 // 13 // Copyright (C) 2000-2008, Intel Corporation, all rights reserved. 14 // Copyright (C) 2009, Willow Garage Inc., all rights reserved. 15 // Third party copyrights are property of their respective owners. 16 // 17 // Redistribution and use in source and binary forms, with or without modification, 18 // are permitted provided that the following conditions are met: 19 // 20 // * Redistribution's of source code must retain the above copyright notice, 21 // this list of conditions and the following disclaimer. 22 // 23 // * Redistribution's in binary form must reproduce the above copyright notice, 24 // this list of conditions and the following disclaimer in the documentation 25 // and/or other materials provided with the distribution. 26 // 27 // * The name of the copyright holders may not be used to endorse or promote products 28 // derived from this software without specific prior written permission. 29 // 30 // This software is provided by the copyright holders and contributors "as is" and 31 // any express or implied warranties, including, but not limited to, the implied 32 // warranties of merchantability and fitness for a particular purpose are disclaimed. 33 // In no event shall the Intel Corporation or contributors be liable for any direct, 34 // indirect, incidental, special, exemplary, or consequential damages 35 // (including, but not limited to, procurement of substitute goods or services; 36 // loss of use, data, or profits; or business interruption) however caused 37 // and on any theory of liability, whether in contract, strict liability, 38 // or tort (including negligence or otherwise) arising in any way out of 39 // the use of this software, even if advised of the possibility of such damage. 40 // 41 //M*/ 42 43 #include "precomp.hpp" 44 #include <limits> 45 46 using namespace cv; 47 48 template<typename _Tp> static int solveQuadratic(_Tp a, _Tp b, _Tp c, _Tp& x1, _Tp& x2) 49 { 50 if( a == 0 ) 51 { 52 if( b == 0 ) 53 { 54 x1 = x2 = 0; 55 return c == 0; 56 } 57 x1 = x2 = -c/b; 58 return 1; 59 } 60 61 _Tp d = b*b - 4*a*c; 62 if( d < 0 ) 63 { 64 x1 = x2 = 0; 65 return 0; 66 } 67 if( d > 0 ) 68 { 69 d = std::sqrt(d); 70 double s = 1/(2*a); 71 x1 = (-b - d)*s; 72 x2 = (-b + d)*s; 73 if( x1 > x2 ) 74 std::swap(x1, x2); 75 return 2; 76 } 77 x1 = x2 = -b/(2*a); 78 return 1; 79 } 80 81 //for android ndk 82 #undef _S 83 static inline Point2f applyHomography( const Mat_<double>& H, const Point2f& pt ) 84 { 85 double z = H(2,0)*pt.x + H(2,1)*pt.y + H(2,2); 86 if( z ) 87 { 88 double w = 1./z; 89 return Point2f( (float)((H(0,0)*pt.x + H(0,1)*pt.y + H(0,2))*w), (float)((H(1,0)*pt.x + H(1,1)*pt.y + H(1,2))*w) ); 90 } 91 return Point2f( std::numeric_limits<float>::max(), std::numeric_limits<float>::max() ); 92 } 93 94 static inline void linearizeHomographyAt( const Mat_<double>& H, const Point2f& pt, Mat_<double>& A ) 95 { 96 A.create(2,2); 97 double p1 = H(0,0)*pt.x + H(0,1)*pt.y + H(0,2), 98 p2 = H(1,0)*pt.x + H(1,1)*pt.y + H(1,2), 99 p3 = H(2,0)*pt.x + H(2,1)*pt.y + H(2,2), 100 p3_2 = p3*p3; 101 if( p3 ) 102 { 103 A(0,0) = H(0,0)/p3 - p1*H(2,0)/p3_2; // fxdx 104 A(0,1) = H(0,1)/p3 - p1*H(2,1)/p3_2; // fxdy 105 106 A(1,0) = H(1,0)/p3 - p2*H(2,0)/p3_2; // fydx 107 A(1,1) = H(1,1)/p3 - p2*H(2,1)/p3_2; // fydx 108 } 109 else 110 A.setTo(Scalar::all(std::numeric_limits<double>::max())); 111 } 112 113 class EllipticKeyPoint 114 { 115 public: 116 EllipticKeyPoint(); 117 EllipticKeyPoint( const Point2f& _center, const Scalar& _ellipse ); 118 119 static void convert( const std::vector<KeyPoint>& src, std::vector<EllipticKeyPoint>& dst ); 120 static void convert( const std::vector<EllipticKeyPoint>& src, std::vector<KeyPoint>& dst ); 121 122 static Mat_<double> getSecondMomentsMatrix( const Scalar& _ellipse ); 123 Mat_<double> getSecondMomentsMatrix() const; 124 125 void calcProjection( const Mat_<double>& H, EllipticKeyPoint& projection ) const; 126 static void calcProjection( const std::vector<EllipticKeyPoint>& src, const Mat_<double>& H, std::vector<EllipticKeyPoint>& dst ); 127 128 Point2f center; 129 Scalar ellipse; // 3 elements a, b, c: ax^2+2bxy+cy^2=1 130 Size_<float> axes; // half length of ellipse axes 131 Size_<float> boundingBox; // half sizes of bounding box which sides are parallel to the coordinate axes 132 }; 133 134 EllipticKeyPoint::EllipticKeyPoint() 135 { 136 *this = EllipticKeyPoint(Point2f(0,0), Scalar(1, 0, 1) ); 137 } 138 139 EllipticKeyPoint::EllipticKeyPoint( const Point2f& _center, const Scalar& _ellipse ) 140 { 141 center = _center; 142 ellipse = _ellipse; 143 144 double a = ellipse[0], b = ellipse[1], c = ellipse[2]; 145 double ac_b2 = a*c - b*b; 146 double x1, x2; 147 solveQuadratic(1., -(a+c), ac_b2, x1, x2); 148 axes.width = (float)(1/sqrt(x1)); 149 axes.height = (float)(1/sqrt(x2)); 150 151 boundingBox.width = (float)sqrt(ellipse[2]/ac_b2); 152 boundingBox.height = (float)sqrt(ellipse[0]/ac_b2); 153 } 154 155 Mat_<double> EllipticKeyPoint::getSecondMomentsMatrix( const Scalar& _ellipse ) 156 { 157 Mat_<double> M(2, 2); 158 M(0,0) = _ellipse[0]; 159 M(1,0) = M(0,1) = _ellipse[1]; 160 M(1,1) = _ellipse[2]; 161 return M; 162 } 163 164 Mat_<double> EllipticKeyPoint::getSecondMomentsMatrix() const 165 { 166 return getSecondMomentsMatrix(ellipse); 167 } 168 169 void EllipticKeyPoint::calcProjection( const Mat_<double>& H, EllipticKeyPoint& projection ) const 170 { 171 Point2f dstCenter = applyHomography(H, center); 172 173 Mat_<double> invM; invert(getSecondMomentsMatrix(), invM); 174 Mat_<double> Aff; linearizeHomographyAt(H, center, Aff); 175 Mat_<double> dstM; invert(Aff*invM*Aff.t(), dstM); 176 177 projection = EllipticKeyPoint( dstCenter, Scalar(dstM(0,0), dstM(0,1), dstM(1,1)) ); 178 } 179 180 void EllipticKeyPoint::convert( const std::vector<KeyPoint>& src, std::vector<EllipticKeyPoint>& dst ) 181 { 182 if( !src.empty() ) 183 { 184 dst.resize(src.size()); 185 for( size_t i = 0; i < src.size(); i++ ) 186 { 187 float rad = src[i].size/2; 188 CV_Assert( rad ); 189 float fac = 1.f/(rad*rad); 190 dst[i] = EllipticKeyPoint( src[i].pt, Scalar(fac, 0, fac) ); 191 } 192 } 193 } 194 195 void EllipticKeyPoint::convert( const std::vector<EllipticKeyPoint>& src, std::vector<KeyPoint>& dst ) 196 { 197 if( !src.empty() ) 198 { 199 dst.resize(src.size()); 200 for( size_t i = 0; i < src.size(); i++ ) 201 { 202 Size_<float> axes = src[i].axes; 203 float rad = sqrt(axes.height*axes.width); 204 dst[i] = KeyPoint(src[i].center, 2*rad ); 205 } 206 } 207 } 208 209 void EllipticKeyPoint::calcProjection( const std::vector<EllipticKeyPoint>& src, const Mat_<double>& H, std::vector<EllipticKeyPoint>& dst ) 210 { 211 if( !src.empty() ) 212 { 213 CV_Assert( !H.empty() && H.cols == 3 && H.rows == 3); 214 dst.resize(src.size()); 215 std::vector<EllipticKeyPoint>::const_iterator srcIt = src.begin(); 216 std::vector<EllipticKeyPoint>::iterator dstIt = dst.begin(); 217 for( ; srcIt != src.end(); ++srcIt, ++dstIt ) 218 srcIt->calcProjection(H, *dstIt); 219 } 220 } 221 222 static void filterEllipticKeyPointsByImageSize( std::vector<EllipticKeyPoint>& keypoints, const Size& imgSize ) 223 { 224 if( !keypoints.empty() ) 225 { 226 std::vector<EllipticKeyPoint> filtered; 227 filtered.reserve(keypoints.size()); 228 std::vector<EllipticKeyPoint>::const_iterator it = keypoints.begin(); 229 for( int i = 0; it != keypoints.end(); ++it, i++ ) 230 { 231 if( it->center.x + it->boundingBox.width < imgSize.width && 232 it->center.x - it->boundingBox.width > 0 && 233 it->center.y + it->boundingBox.height < imgSize.height && 234 it->center.y - it->boundingBox.height > 0 ) 235 filtered.push_back(*it); 236 } 237 keypoints.assign(filtered.begin(), filtered.end()); 238 } 239 } 240 241 struct IntersectAreaCounter 242 { 243 IntersectAreaCounter( float _dr, int _minx, 244 int _miny, int _maxy, 245 const Point2f& _diff, 246 const Scalar& _ellipse1, const Scalar& _ellipse2 ) : 247 dr(_dr), bua(0), bna(0), minx(_minx), miny(_miny), maxy(_maxy), 248 diff(_diff), ellipse1(_ellipse1), ellipse2(_ellipse2) {} 249 IntersectAreaCounter( const IntersectAreaCounter& counter, Split ) 250 { 251 *this = counter; 252 bua = 0; 253 bna = 0; 254 } 255 256 void operator()( const BlockedRange& range ) 257 { 258 CV_Assert( miny < maxy ); 259 CV_Assert( dr > FLT_EPSILON ); 260 261 int temp_bua = bua, temp_bna = bna; 262 for( int i = range.begin(); i != range.end(); i++ ) 263 { 264 float rx1 = minx + i*dr; 265 float rx2 = rx1 - diff.x; 266 for( float ry1 = (float)miny; ry1 <= (float)maxy; ry1 += dr ) 267 { 268 float ry2 = ry1 - diff.y; 269 //compute the distance from the ellipse center 270 float e1 = (float)(ellipse1[0]*rx1*rx1 + 2*ellipse1[1]*rx1*ry1 + ellipse1[2]*ry1*ry1); 271 float e2 = (float)(ellipse2[0]*rx2*rx2 + 2*ellipse2[1]*rx2*ry2 + ellipse2[2]*ry2*ry2); 272 //compute the area 273 if( e1<1 && e2<1 ) temp_bna++; 274 if( e1<1 || e2<1 ) temp_bua++; 275 } 276 } 277 bua = temp_bua; 278 bna = temp_bna; 279 } 280 281 void join( IntersectAreaCounter& ac ) 282 { 283 bua += ac.bua; 284 bna += ac.bna; 285 } 286 287 float dr; 288 int bua, bna; 289 290 int minx; 291 int miny, maxy; 292 293 Point2f diff; 294 Scalar ellipse1, ellipse2; 295 296 }; 297 298 struct SIdx 299 { 300 SIdx() : S(-1), i1(-1), i2(-1) {} 301 SIdx(float _S, int _i1, int _i2) : S(_S), i1(_i1), i2(_i2) {} 302 float S; 303 int i1; 304 int i2; 305 306 bool operator<(const SIdx& v) const { return S > v.S; } 307 308 struct UsedFinder 309 { 310 UsedFinder(const SIdx& _used) : used(_used) {} 311 const SIdx& used; 312 bool operator()(const SIdx& v) const { return (v.i1 == used.i1 || v.i2 == used.i2); } 313 UsedFinder& operator=(const UsedFinder&); 314 }; 315 }; 316 317 static void computeOneToOneMatchedOverlaps( const std::vector<EllipticKeyPoint>& keypoints1, const std::vector<EllipticKeyPoint>& keypoints2t, 318 bool commonPart, std::vector<SIdx>& overlaps, float minOverlap ) 319 { 320 CV_Assert( minOverlap >= 0.f ); 321 overlaps.clear(); 322 if( keypoints1.empty() || keypoints2t.empty() ) 323 return; 324 325 overlaps.clear(); 326 overlaps.reserve(cvRound(keypoints1.size() * keypoints2t.size() * 0.01)); 327 328 for( size_t i1 = 0; i1 < keypoints1.size(); i1++ ) 329 { 330 EllipticKeyPoint kp1 = keypoints1[i1]; 331 float maxDist = sqrt(kp1.axes.width*kp1.axes.height), 332 fac = 30.f/maxDist; 333 if( !commonPart ) 334 fac=3; 335 336 maxDist = maxDist*4; 337 fac = 1.f/(fac*fac); 338 339 EllipticKeyPoint keypoint1a = EllipticKeyPoint( kp1.center, Scalar(fac*kp1.ellipse[0], fac*kp1.ellipse[1], fac*kp1.ellipse[2]) ); 340 341 for( size_t i2 = 0; i2 < keypoints2t.size(); i2++ ) 342 { 343 EllipticKeyPoint kp2 = keypoints2t[i2]; 344 Point2f diff = kp2.center - kp1.center; 345 346 if( norm(diff) < maxDist ) 347 { 348 EllipticKeyPoint keypoint2a = EllipticKeyPoint( kp2.center, Scalar(fac*kp2.ellipse[0], fac*kp2.ellipse[1], fac*kp2.ellipse[2]) ); 349 //find the largest eigenvalue 350 int maxx = (int)ceil(( keypoint1a.boundingBox.width > (diff.x+keypoint2a.boundingBox.width)) ? 351 keypoint1a.boundingBox.width : (diff.x+keypoint2a.boundingBox.width)); 352 int minx = (int)floor((-keypoint1a.boundingBox.width < (diff.x-keypoint2a.boundingBox.width)) ? 353 -keypoint1a.boundingBox.width : (diff.x-keypoint2a.boundingBox.width)); 354 355 int maxy = (int)ceil(( keypoint1a.boundingBox.height > (diff.y+keypoint2a.boundingBox.height)) ? 356 keypoint1a.boundingBox.height : (diff.y+keypoint2a.boundingBox.height)); 357 int miny = (int)floor((-keypoint1a.boundingBox.height < (diff.y-keypoint2a.boundingBox.height)) ? 358 -keypoint1a.boundingBox.height : (diff.y-keypoint2a.boundingBox.height)); 359 int mina = (maxx-minx) < (maxy-miny) ? (maxx-minx) : (maxy-miny) ; 360 361 //compute the area 362 float dr = (float)mina/50.f; 363 int N = (int)floor((float)(maxx - minx) / dr); 364 IntersectAreaCounter ac( dr, minx, miny, maxy, diff, keypoint1a.ellipse, keypoint2a.ellipse ); 365 parallel_reduce( BlockedRange(0, N+1), ac ); 366 if( ac.bna > 0 ) 367 { 368 float ov = (float)ac.bna / (float)ac.bua; 369 if( ov >= minOverlap ) 370 overlaps.push_back(SIdx(ov, (int)i1, (int)i2)); 371 } 372 } 373 } 374 } 375 376 std::sort( overlaps.begin(), overlaps.end() ); 377 378 typedef std::vector<SIdx>::iterator It; 379 380 It pos = overlaps.begin(); 381 It end = overlaps.end(); 382 383 while(pos != end) 384 { 385 It prev = pos++; 386 end = std::remove_if(pos, end, SIdx::UsedFinder(*prev)); 387 } 388 overlaps.erase(pos, overlaps.end()); 389 } 390 391 static void calculateRepeatability( const Mat& img1, const Mat& img2, const Mat& H1to2, 392 const std::vector<KeyPoint>& _keypoints1, const std::vector<KeyPoint>& _keypoints2, 393 float& repeatability, int& correspondencesCount, 394 Mat* thresholdedOverlapMask=0 ) 395 { 396 std::vector<EllipticKeyPoint> keypoints1, keypoints2, keypoints1t, keypoints2t; 397 EllipticKeyPoint::convert( _keypoints1, keypoints1 ); 398 EllipticKeyPoint::convert( _keypoints2, keypoints2 ); 399 400 // calculate projections of key points 401 EllipticKeyPoint::calcProjection( keypoints1, H1to2, keypoints1t ); 402 Mat H2to1; invert(H1to2, H2to1); 403 EllipticKeyPoint::calcProjection( keypoints2, H2to1, keypoints2t ); 404 405 float overlapThreshold; 406 bool ifEvaluateDetectors = thresholdedOverlapMask == 0; 407 if( ifEvaluateDetectors ) 408 { 409 overlapThreshold = 1.f - 0.4f; 410 411 // remove key points from outside of the common image part 412 Size sz1 = img1.size(), sz2 = img2.size(); 413 filterEllipticKeyPointsByImageSize( keypoints1, sz1 ); 414 filterEllipticKeyPointsByImageSize( keypoints1t, sz2 ); 415 filterEllipticKeyPointsByImageSize( keypoints2, sz2 ); 416 filterEllipticKeyPointsByImageSize( keypoints2t, sz1 ); 417 } 418 else 419 { 420 overlapThreshold = 1.f - 0.5f; 421 422 thresholdedOverlapMask->create( (int)keypoints1.size(), (int)keypoints2t.size(), CV_8UC1 ); 423 thresholdedOverlapMask->setTo( Scalar::all(0) ); 424 } 425 size_t size1 = keypoints1.size(), size2 = keypoints2t.size(); 426 size_t minCount = MIN( size1, size2 ); 427 428 // calculate overlap errors 429 std::vector<SIdx> overlaps; 430 computeOneToOneMatchedOverlaps( keypoints1, keypoints2t, ifEvaluateDetectors, overlaps, overlapThreshold/*min overlap*/ ); 431 432 correspondencesCount = -1; 433 repeatability = -1.f; 434 if( overlaps.empty() ) 435 return; 436 437 if( ifEvaluateDetectors ) 438 { 439 // regions one-to-one matching 440 correspondencesCount = (int)overlaps.size(); 441 repeatability = minCount ? (float)correspondencesCount / minCount : -1; 442 } 443 else 444 { 445 for( size_t i = 0; i < overlaps.size(); i++ ) 446 { 447 int y = overlaps[i].i1; 448 int x = overlaps[i].i2; 449 thresholdedOverlapMask->at<uchar>(y,x) = 1; 450 } 451 } 452 } 453 454 void cv::evaluateFeatureDetector( const Mat& img1, const Mat& img2, const Mat& H1to2, 455 std::vector<KeyPoint>* _keypoints1, std::vector<KeyPoint>* _keypoints2, 456 float& repeatability, int& correspCount, 457 const Ptr<FeatureDetector>& _fdetector ) 458 { 459 Ptr<FeatureDetector> fdetector(_fdetector); 460 std::vector<KeyPoint> *keypoints1, *keypoints2, buf1, buf2; 461 keypoints1 = _keypoints1 != 0 ? _keypoints1 : &buf1; 462 keypoints2 = _keypoints2 != 0 ? _keypoints2 : &buf2; 463 464 if( (keypoints1->empty() || keypoints2->empty()) && !fdetector ) 465 CV_Error( Error::StsBadArg, "fdetector must not be empty when keypoints1 or keypoints2 is empty" ); 466 467 if( keypoints1->empty() ) 468 fdetector->detect( img1, *keypoints1 ); 469 if( keypoints2->empty() ) 470 fdetector->detect( img2, *keypoints2 ); 471 472 calculateRepeatability( img1, img2, H1to2, *keypoints1, *keypoints2, repeatability, correspCount ); 473 } 474 475 struct DMatchForEvaluation : public DMatch 476 { 477 uchar isCorrect; 478 DMatchForEvaluation( const DMatch &dm ) : DMatch( dm ) {} 479 }; 480 481 static inline float recall( int correctMatchCount, int correspondenceCount ) 482 { 483 return correspondenceCount ? (float)correctMatchCount / (float)correspondenceCount : -1; 484 } 485 486 static inline float precision( int correctMatchCount, int falseMatchCount ) 487 { 488 return correctMatchCount + falseMatchCount ? (float)correctMatchCount / (float)(correctMatchCount + falseMatchCount) : -1; 489 } 490 491 void cv::computeRecallPrecisionCurve( const std::vector<std::vector<DMatch> >& matches1to2, 492 const std::vector<std::vector<uchar> >& correctMatches1to2Mask, 493 std::vector<Point2f>& recallPrecisionCurve ) 494 { 495 CV_Assert( matches1to2.size() == correctMatches1to2Mask.size() ); 496 497 std::vector<DMatchForEvaluation> allMatches; 498 int correspondenceCount = 0; 499 for( size_t i = 0; i < matches1to2.size(); i++ ) 500 { 501 for( size_t j = 0; j < matches1to2[i].size(); j++ ) 502 { 503 DMatchForEvaluation match = matches1to2[i][j]; 504 match.isCorrect = correctMatches1to2Mask[i][j] ; 505 allMatches.push_back( match ); 506 correspondenceCount += match.isCorrect != 0 ? 1 : 0; 507 } 508 } 509 510 std::sort( allMatches.begin(), allMatches.end() ); 511 512 int correctMatchCount = 0, falseMatchCount = 0; 513 recallPrecisionCurve.resize( allMatches.size() ); 514 for( size_t i = 0; i < allMatches.size(); i++ ) 515 { 516 if( allMatches[i].isCorrect ) 517 correctMatchCount++; 518 else 519 falseMatchCount++; 520 521 float r = recall( correctMatchCount, correspondenceCount ); 522 float p = precision( correctMatchCount, falseMatchCount ); 523 recallPrecisionCurve[i] = Point2f(1-p, r); 524 } 525 } 526 527 float cv::getRecall( const std::vector<Point2f>& recallPrecisionCurve, float l_precision ) 528 { 529 int nearestPointIndex = getNearestPoint( recallPrecisionCurve, l_precision ); 530 531 float recall = -1.f; 532 533 if( nearestPointIndex >= 0 ) 534 recall = recallPrecisionCurve[nearestPointIndex].y; 535 536 return recall; 537 } 538 539 int cv::getNearestPoint( const std::vector<Point2f>& recallPrecisionCurve, float l_precision ) 540 { 541 int nearestPointIndex = -1; 542 543 if( l_precision >= 0 && l_precision <= 1 ) 544 { 545 float minDiff = FLT_MAX; 546 for( size_t i = 0; i < recallPrecisionCurve.size(); i++ ) 547 { 548 float curDiff = std::fabs(l_precision - recallPrecisionCurve[i].x); 549 if( curDiff <= minDiff ) 550 { 551 nearestPointIndex = (int)i; 552 minDiff = curDiff; 553 } 554 } 555 } 556 557 return nearestPointIndex; 558 } 559
