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 /* 44 * Implementation of the paper Shape Matching and Object Recognition Using Shape Contexts 45 * Belongie et al., 2002 by Juan Manuel Perez for GSoC 2013. 46 */ 47 48 #include "precomp.hpp" 49 #include "opencv2/core.hpp" 50 #include "scd_def.hpp" 51 #include <limits> 52 53 namespace cv 54 { 55 class ShapeContextDistanceExtractorImpl : public ShapeContextDistanceExtractor 56 { 57 public: 58 /* Constructors */ 59 ShapeContextDistanceExtractorImpl(int _nAngularBins, int _nRadialBins, float _innerRadius, float _outerRadius, int _iterations, 60 const Ptr<HistogramCostExtractor> &_comparer, const Ptr<ShapeTransformer> &_transformer) 61 { 62 nAngularBins=_nAngularBins; 63 nRadialBins=_nRadialBins; 64 innerRadius=_innerRadius; 65 outerRadius=_outerRadius; 66 rotationInvariant=false; 67 comparer=_comparer; 68 iterations=_iterations; 69 transformer=_transformer; 70 bendingEnergyWeight=0.3f; 71 imageAppearanceWeight=0.0f; 72 shapeContextWeight=1.0f; 73 sigma=10.0f; 74 name_ = "ShapeDistanceExtractor.SCD"; 75 } 76 77 /* Destructor */ 78 ~ShapeContextDistanceExtractorImpl() 79 { 80 } 81 82 //! the main operator 83 virtual float computeDistance(InputArray contour1, InputArray contour2); 84 85 //! Setters/Getters 86 virtual void setAngularBins(int _nAngularBins){CV_Assert(_nAngularBins>0); nAngularBins=_nAngularBins;} 87 virtual int getAngularBins() const {return nAngularBins;} 88 89 virtual void setRadialBins(int _nRadialBins){CV_Assert(_nRadialBins>0); nRadialBins=_nRadialBins;} 90 virtual int getRadialBins() const {return nRadialBins;} 91 92 virtual void setInnerRadius(float _innerRadius) {CV_Assert(_innerRadius>0); innerRadius=_innerRadius;} 93 virtual float getInnerRadius() const {return innerRadius;} 94 95 virtual void setOuterRadius(float _outerRadius) {CV_Assert(_outerRadius>0); outerRadius=_outerRadius;} 96 virtual float getOuterRadius() const {return outerRadius;} 97 98 virtual void setRotationInvariant(bool _rotationInvariant) {rotationInvariant=_rotationInvariant;} 99 virtual bool getRotationInvariant() const {return rotationInvariant;} 100 101 virtual void setCostExtractor(Ptr<HistogramCostExtractor> _comparer) { comparer = _comparer; } 102 virtual Ptr<HistogramCostExtractor> getCostExtractor() const { return comparer; } 103 104 virtual void setShapeContextWeight(float _shapeContextWeight) {shapeContextWeight=_shapeContextWeight;} 105 virtual float getShapeContextWeight() const {return shapeContextWeight;} 106 107 virtual void setImageAppearanceWeight(float _imageAppearanceWeight) {imageAppearanceWeight=_imageAppearanceWeight;} 108 virtual float getImageAppearanceWeight() const {return imageAppearanceWeight;} 109 110 virtual void setBendingEnergyWeight(float _bendingEnergyWeight) {bendingEnergyWeight=_bendingEnergyWeight;} 111 virtual float getBendingEnergyWeight() const {return bendingEnergyWeight;} 112 113 virtual void setStdDev(float _sigma) {sigma=_sigma;} 114 virtual float getStdDev() const {return sigma;} 115 116 virtual void setImages(InputArray _image1, InputArray _image2) 117 { 118 Mat image1_=_image1.getMat(), image2_=_image2.getMat(); 119 CV_Assert((image1_.depth()==0) && (image2_.depth()==0)); 120 image1=image1_; 121 image2=image2_; 122 } 123 124 virtual void getImages(OutputArray _image1, OutputArray _image2) const 125 { 126 CV_Assert((!image1.empty()) && (!image2.empty())); 127 _image1.create(image1.size(), image1.type()); 128 _image2.create(image2.size(), image2.type()); 129 _image1.getMat()=image1; 130 _image2.getMat()=image2; 131 } 132 133 virtual void setIterations(int _iterations) {CV_Assert(_iterations>0); iterations=_iterations;} 134 virtual int getIterations() const {return iterations;} 135 136 virtual void setTransformAlgorithm(Ptr<ShapeTransformer> _transformer) {transformer=_transformer;} 137 virtual Ptr<ShapeTransformer> getTransformAlgorithm() const {return transformer;} 138 139 //! write/read 140 virtual void write(FileStorage& fs) const 141 { 142 fs << "name" << name_ 143 << "nRads" << nRadialBins 144 << "nAngs" << nAngularBins 145 << "iters" << iterations 146 << "img_1" << image1 147 << "img_2" << image2 148 << "beWei" << bendingEnergyWeight 149 << "scWei" << shapeContextWeight 150 << "iaWei" << imageAppearanceWeight 151 << "costF" << costFlag 152 << "rotIn" << rotationInvariant 153 << "sigma" << sigma; 154 } 155 156 virtual void read(const FileNode& fn) 157 { 158 CV_Assert( (String)fn["name"] == name_ ); 159 nRadialBins = (int)fn["nRads"]; 160 nAngularBins = (int)fn["nAngs"]; 161 iterations = (int)fn["iters"]; 162 bendingEnergyWeight = (float)fn["beWei"]; 163 shapeContextWeight = (float)fn["scWei"]; 164 imageAppearanceWeight = (float)fn["iaWei"]; 165 costFlag = (int)fn["costF"]; 166 sigma = (float)fn["sigma"]; 167 } 168 169 protected: 170 int nAngularBins; 171 int nRadialBins; 172 float innerRadius; 173 float outerRadius; 174 bool rotationInvariant; 175 int costFlag; 176 int iterations; 177 Ptr<ShapeTransformer> transformer; 178 Ptr<HistogramCostExtractor> comparer; 179 Mat image1; 180 Mat image2; 181 float bendingEnergyWeight; 182 float imageAppearanceWeight; 183 float shapeContextWeight; 184 float sigma; 185 String name_; 186 }; 187 188 float ShapeContextDistanceExtractorImpl::computeDistance(InputArray contour1, InputArray contour2) 189 { 190 // Checking // 191 Mat sset1=contour1.getMat(), sset2=contour2.getMat(), set1, set2; 192 if (set1.type() != CV_32F) 193 sset1.convertTo(set1, CV_32F); 194 else 195 sset1.copyTo(set1); 196 197 if (set2.type() != CV_32F) 198 sset2.convertTo(set2, CV_32F); 199 else 200 sset2.copyTo(set2); 201 202 CV_Assert((set1.channels()==2) && (set1.cols>0)); 203 CV_Assert((set2.channels()==2) && (set2.cols>0)); 204 if (imageAppearanceWeight!=0) 205 { 206 CV_Assert((!image1.empty()) && (!image2.empty())); 207 } 208 209 // Initializing Extractor, Descriptor structures and Matcher // 210 SCD set1SCE(nAngularBins, nRadialBins, innerRadius, outerRadius, rotationInvariant); 211 Mat set1SCD; 212 SCD set2SCE(nAngularBins, nRadialBins, innerRadius, outerRadius, rotationInvariant); 213 Mat set2SCD; 214 SCDMatcher matcher; 215 std::vector<DMatch> matches; 216 217 // Distance components (The output is a linear combination of these 3) // 218 float sDistance=0, bEnergy=0, iAppearance=0; 219 float beta; 220 221 // Initializing some variables // 222 std::vector<int> inliers1, inliers2; 223 224 Ptr<ThinPlateSplineShapeTransformer> transDown = transformer.dynamicCast<ThinPlateSplineShapeTransformer>(); 225 226 Mat warpedImage; 227 int ii, jj, pt; 228 229 for (ii=0; ii<iterations; ii++) 230 { 231 // Extract SCD descriptor in the set1 // 232 set1SCE.extractSCD(set1, set1SCD, inliers1); 233 234 // Extract SCD descriptor of the set2 (TARGET) // 235 set2SCE.extractSCD(set2, set2SCD, inliers2, set1SCE.getMeanDistance()); 236 237 // regularization parameter with annealing rate annRate // 238 beta=set1SCE.getMeanDistance(); 239 beta *= beta; 240 241 // match // 242 matcher.matchDescriptors(set1SCD, set2SCD, matches, comparer, inliers1, inliers2); 243 244 // apply TPS transform // 245 if ( !transDown.empty() ) 246 transDown->setRegularizationParameter(beta); 247 transformer->estimateTransformation(set1, set2, matches); 248 bEnergy += transformer->applyTransformation(set1, set1); 249 250 // Image appearance // 251 if (imageAppearanceWeight!=0) 252 { 253 // Have to accumulate the transformation along all the iterations 254 if (ii==0) 255 { 256 if ( !transDown.empty() ) 257 { 258 image2.copyTo(warpedImage); 259 } 260 else 261 { 262 image1.copyTo(warpedImage); 263 } 264 } 265 transformer->warpImage(warpedImage, warpedImage); 266 } 267 } 268 269 Mat gaussWindow, diffIm; 270 if (imageAppearanceWeight!=0) 271 { 272 // compute appearance cost 273 if ( !transDown.empty() ) 274 { 275 resize(warpedImage, warpedImage, image1.size()); 276 Mat temp=(warpedImage-image1); 277 multiply(temp, temp, diffIm); 278 } 279 else 280 { 281 resize(warpedImage, warpedImage, image2.size()); 282 Mat temp=(warpedImage-image2); 283 multiply(temp, temp, diffIm); 284 } 285 gaussWindow = Mat::zeros(warpedImage.rows, warpedImage.cols, CV_32F); 286 for (pt=0; pt<sset1.cols; pt++) 287 { 288 Point2f p = sset1.at<Point2f>(0,pt); 289 for (ii=0; ii<diffIm.rows; ii++) 290 { 291 for (jj=0; jj<diffIm.cols; jj++) 292 { 293 float val = float(std::exp( -float( (p.x-jj)*(p.x-jj) + (p.y-ii)*(p.y-ii) )/(2*sigma*sigma) ) / (sigma*sigma*2*CV_PI)); 294 gaussWindow.at<float>(ii,jj) += val; 295 } 296 } 297 } 298 299 Mat appIm(diffIm.rows, diffIm.cols, CV_32F); 300 for (ii=0; ii<diffIm.rows; ii++) 301 { 302 for (jj=0; jj<diffIm.cols; jj++) 303 { 304 float elema=float( diffIm.at<uchar>(ii,jj) )/255; 305 float elemb=gaussWindow.at<float>(ii,jj); 306 appIm.at<float>(ii,jj) = elema*elemb; 307 } 308 } 309 iAppearance = float(cv::sum(appIm)[0]/sset1.cols); 310 } 311 sDistance = matcher.getMatchingCost(); 312 313 return (sDistance*shapeContextWeight+bEnergy*bendingEnergyWeight+iAppearance*imageAppearanceWeight); 314 } 315 316 Ptr <ShapeContextDistanceExtractor> createShapeContextDistanceExtractor(int nAngularBins, int nRadialBins, float innerRadius, float outerRadius, int iterations, 317 const Ptr<HistogramCostExtractor> &comparer, const Ptr<ShapeTransformer> &transformer) 318 { 319 return Ptr <ShapeContextDistanceExtractor> ( new ShapeContextDistanceExtractorImpl(nAngularBins, nRadialBins, innerRadius, 320 outerRadius, iterations, comparer, transformer) ); 321 } 322 323 //! SCD 324 void SCD::extractSCD(cv::Mat &contour, cv::Mat &descriptors, const std::vector<int> &queryInliers, const float _meanDistance) 325 { 326 cv::Mat contourMat = contour; 327 cv::Mat disMatrix = cv::Mat::zeros(contourMat.cols, contourMat.cols, CV_32F); 328 cv::Mat angleMatrix = cv::Mat::zeros(contourMat.cols, contourMat.cols, CV_32F); 329 330 std::vector<double> logspaces, angspaces; 331 logarithmicSpaces(logspaces); 332 angularSpaces(angspaces); 333 buildNormalizedDistanceMatrix(contourMat, disMatrix, queryInliers, _meanDistance); 334 buildAngleMatrix(contourMat, angleMatrix); 335 336 // Now, build the descriptor matrix (each row is a point) // 337 descriptors = cv::Mat::zeros(contourMat.cols, descriptorSize(), CV_32F); 338 339 for (int ptidx=0; ptidx<contourMat.cols; ptidx++) 340 { 341 for (int cmp=0; cmp<contourMat.cols; cmp++) 342 { 343 if (ptidx==cmp) continue; 344 if ((int)queryInliers.size()>0) 345 { 346 if (queryInliers[ptidx]==0 || queryInliers[cmp]==0) continue; //avoid outliers 347 } 348 349 int angidx=-1, radidx=-1; 350 for (int i=0; i<nRadialBins; i++) 351 { 352 if (disMatrix.at<float>(ptidx, cmp)<logspaces[i]) 353 { 354 radidx=i; 355 break; 356 } 357 } 358 for (int i=0; i<nAngularBins; i++) 359 { 360 if (angleMatrix.at<float>(ptidx, cmp)<angspaces[i]) 361 { 362 angidx=i; 363 break; 364 } 365 } 366 if (angidx!=-1 && radidx!=-1) 367 { 368 int idx = angidx+radidx*nAngularBins; 369 descriptors.at<float>(ptidx, idx)++; 370 } 371 } 372 } 373 } 374 375 void SCD::logarithmicSpaces(std::vector<double> &vecSpaces) const 376 { 377 double logmin=log10(innerRadius); 378 double logmax=log10(outerRadius); 379 double delta=(logmax-logmin)/(nRadialBins-1); 380 double accdelta=0; 381 382 for (int i=0; i<nRadialBins; i++) 383 { 384 double val = std::pow(10,logmin+accdelta); 385 vecSpaces.push_back(val); 386 accdelta += delta; 387 } 388 } 389 390 void SCD::angularSpaces(std::vector<double> &vecSpaces) const 391 { 392 double delta=2*CV_PI/nAngularBins; 393 double val=0; 394 395 for (int i=0; i<nAngularBins; i++) 396 { 397 val += delta; 398 vecSpaces.push_back(val); 399 } 400 } 401 402 void SCD::buildNormalizedDistanceMatrix(cv::Mat &contour, cv::Mat &disMatrix, const std::vector<int> &queryInliers, const float _meanDistance) 403 { 404 cv::Mat contourMat = contour; 405 cv::Mat mask(disMatrix.rows, disMatrix.cols, CV_8U); 406 407 for (int i=0; i<contourMat.cols; i++) 408 { 409 for (int j=0; j<contourMat.cols; j++) 410 { 411 disMatrix.at<float>(i,j) = (float)norm( cv::Mat(contourMat.at<cv::Point2f>(0,i)-contourMat.at<cv::Point2f>(0,j)), cv::NORM_L2 ); 412 if (_meanDistance<0) 413 { 414 if (queryInliers.size()>0) 415 { 416 mask.at<char>(i,j)=char(queryInliers[j] && queryInliers[i]); 417 } 418 else 419 { 420 mask.at<char>(i,j)=1; 421 } 422 } 423 } 424 } 425 426 if (_meanDistance<0) 427 { 428 meanDistance=(float)mean(disMatrix, mask)[0]; 429 } 430 else 431 { 432 meanDistance=_meanDistance; 433 } 434 disMatrix/=meanDistance+FLT_EPSILON; 435 } 436 437 void SCD::buildAngleMatrix(cv::Mat &contour, cv::Mat &angleMatrix) const 438 { 439 cv::Mat contourMat = contour; 440 441 // if descriptor is rotationInvariant compute massCenter // 442 cv::Point2f massCenter(0,0); 443 if (rotationInvariant) 444 { 445 for (int i=0; i<contourMat.cols; i++) 446 { 447 massCenter+=contourMat.at<cv::Point2f>(0,i); 448 } 449 massCenter.x=massCenter.x/(float)contourMat.cols; 450 massCenter.y=massCenter.y/(float)contourMat.cols; 451 } 452 453 454 for (int i=0; i<contourMat.cols; i++) 455 { 456 for (int j=0; j<contourMat.cols; j++) 457 { 458 if (i==j) 459 { 460 angleMatrix.at<float>(i,j)=0.0; 461 } 462 else 463 { 464 cv::Point2f dif = contourMat.at<cv::Point2f>(0,i) - contourMat.at<cv::Point2f>(0,j); 465 angleMatrix.at<float>(i,j) = std::atan2(dif.y, dif.x); 466 467 if (rotationInvariant) 468 { 469 cv::Point2f refPt = contourMat.at<cv::Point2f>(0,i) - massCenter; 470 float refAngle = atan2(refPt.y, refPt.x); 471 angleMatrix.at<float>(i,j) -= refAngle; 472 } 473 angleMatrix.at<float>(i,j) = float(fmod(double(angleMatrix.at<float>(i,j)+(double)FLT_EPSILON),2*CV_PI)+CV_PI); 474 } 475 } 476 } 477 } 478 479 //! SCDMatcher 480 void SCDMatcher::matchDescriptors(cv::Mat &descriptors1, cv::Mat &descriptors2, std::vector<cv::DMatch> &matches, 481 cv::Ptr<cv::HistogramCostExtractor> &comparer, std::vector<int> &inliers1, std::vector<int> &inliers2) 482 { 483 matches.clear(); 484 485 // Build the cost Matrix between descriptors // 486 cv::Mat costMat; 487 buildCostMatrix(descriptors1, descriptors2, costMat, comparer); 488 489 // Solve the matching problem using the hungarian method // 490 hungarian(costMat, matches, inliers1, inliers2, descriptors1.rows, descriptors2.rows); 491 } 492 493 void SCDMatcher::buildCostMatrix(const cv::Mat &descriptors1, const cv::Mat &descriptors2, 494 cv::Mat &costMatrix, cv::Ptr<cv::HistogramCostExtractor> &comparer) const 495 { 496 comparer->buildCostMatrix(descriptors1, descriptors2, costMatrix); 497 } 498 499 void SCDMatcher::hungarian(cv::Mat &costMatrix, std::vector<cv::DMatch> &outMatches, std::vector<int> &inliers1, 500 std::vector<int> &inliers2, int sizeScd1, int sizeScd2) 501 { 502 std::vector<int> free(costMatrix.rows, 0), collist(costMatrix.rows, 0); 503 std::vector<int> matches(costMatrix.rows, 0), colsol(costMatrix.rows), rowsol(costMatrix.rows); 504 std::vector<float> d(costMatrix.rows), pred(costMatrix.rows), v(costMatrix.rows); 505 506 const float LOWV = 1e-10f; 507 bool unassignedfound; 508 int i=0, imin=0, numfree=0, prvnumfree=0, f=0, i0=0, k=0, freerow=0; 509 int j=0, j1=0, j2=0, endofpath=0, last=0, low=0, up=0; 510 float min=0, h=0, umin=0, usubmin=0, v2=0; 511 512 // COLUMN REDUCTION // 513 for (j = costMatrix.rows-1; j >= 0; j--) 514 { 515 // find minimum cost over rows. 516 min = costMatrix.at<float>(0,j); 517 imin = 0; 518 for (i = 1; i < costMatrix.rows; i++) 519 if (costMatrix.at<float>(i,j) < min) 520 { 521 min = costMatrix.at<float>(i,j); 522 imin = i; 523 } 524 v[j] = min; 525 526 if (++matches[imin] == 1) 527 { 528 rowsol[imin] = j; 529 colsol[j] = imin; 530 } 531 else 532 { 533 colsol[j]=-1; 534 } 535 } 536 537 // REDUCTION TRANSFER // 538 for (i=0; i<costMatrix.rows; i++) 539 { 540 if (matches[i] == 0) 541 { 542 free[numfree++] = i; 543 } 544 else 545 { 546 if (matches[i] == 1) 547 { 548 j1=rowsol[i]; 549 min=std::numeric_limits<float>::max(); 550 for (j=0; j<costMatrix.rows; j++) 551 { 552 if (j!=j1) 553 { 554 if (costMatrix.at<float>(i,j)-v[j] < min) 555 { 556 min=costMatrix.at<float>(i,j)-v[j]; 557 } 558 } 559 } 560 v[j1] = v[j1]-min; 561 } 562 } 563 } 564 // AUGMENTING ROW REDUCTION // 565 int loopcnt = 0; 566 do 567 { 568 loopcnt++; 569 k=0; 570 prvnumfree=numfree; 571 numfree=0; 572 while (k < prvnumfree) 573 { 574 i=free[k]; 575 k++; 576 umin = costMatrix.at<float>(i,0)-v[0]; 577 j1=0; 578 usubmin = std::numeric_limits<float>::max(); 579 for (j=1; j<costMatrix.rows; j++) 580 { 581 h = costMatrix.at<float>(i,j)-v[j]; 582 if (h < usubmin) 583 { 584 if (h >= umin) 585 { 586 usubmin = h; 587 j2 = j; 588 } 589 else 590 { 591 usubmin = umin; 592 umin = h; 593 j2 = j1; 594 j1 = j; 595 } 596 } 597 } 598 i0 = colsol[j1]; 599 600 if (fabs(umin-usubmin) > LOWV) //if( umin < usubmin ) 601 { 602 v[j1] = v[j1] - (usubmin - umin); 603 } 604 else // minimum and subminimum equal. 605 { 606 if (i0 >= 0) // minimum column j1 is assigned. 607 { 608 j1 = j2; 609 i0 = colsol[j2]; 610 } 611 } 612 // (re-)assign i to j1, possibly de-assigning an i0. 613 rowsol[i]=j1; 614 colsol[j1]=i; 615 616 if (i0 >= 0) 617 { 618 //if( umin < usubmin ) 619 if (fabs(umin-usubmin) > LOWV) 620 { 621 free[--k] = i0; 622 } 623 else 624 { 625 free[numfree++] = i0; 626 } 627 } 628 } 629 }while (loopcnt<2); // repeat once. 630 631 // AUGMENT SOLUTION for each free row // 632 for (f = 0; f<numfree; f++) 633 { 634 freerow = free[f]; // start row of augmenting path. 635 // Dijkstra shortest path algorithm. 636 // runs until unassigned column added to shortest path tree. 637 for (j = 0; j < costMatrix.rows; j++) 638 { 639 d[j] = costMatrix.at<float>(freerow,j) - v[j]; 640 pred[j] = float(freerow); 641 collist[j] = j; // init column list. 642 } 643 644 low=0; // columns in 0..low-1 are ready, now none. 645 up=0; // columns in low..up-1 are to be scanned for current minimum, now none. 646 unassignedfound = false; 647 do 648 { 649 if (up == low) 650 { 651 last=low-1; 652 min = d[collist[up++]]; 653 for (k = up; k < costMatrix.rows; k++) 654 { 655 j = collist[k]; 656 h = d[j]; 657 if (h <= min) 658 { 659 if (h < min) // new minimum. 660 { 661 up = low; // restart list at index low. 662 min = h; 663 } 664 collist[k] = collist[up]; 665 collist[up++] = j; 666 } 667 } 668 for (k=low; k<up; k++) 669 { 670 if (colsol[collist[k]] < 0) 671 { 672 endofpath = collist[k]; 673 unassignedfound = true; 674 break; 675 } 676 } 677 } 678 679 if (!unassignedfound) 680 { 681 // update 'distances' between freerow and all unscanned columns, via next scanned column. 682 j1 = collist[low]; 683 low++; 684 i = colsol[j1]; 685 h = costMatrix.at<float>(i,j1)-v[j1]-min; 686 687 for (k = up; k < costMatrix.rows; k++) 688 { 689 j = collist[k]; 690 v2 = costMatrix.at<float>(i,j) - v[j] - h; 691 if (v2 < d[j]) 692 { 693 pred[j] = float(i); 694 if (v2 == min) 695 { 696 if (colsol[j] < 0) 697 { 698 // if unassigned, shortest augmenting path is complete. 699 endofpath = j; 700 unassignedfound = true; 701 break; 702 } 703 else 704 { 705 collist[k] = collist[up]; 706 collist[up++] = j; 707 } 708 } 709 d[j] = v2; 710 } 711 } 712 } 713 }while (!unassignedfound); 714 715 // update column prices. 716 for (k = 0; k <= last; k++) 717 { 718 j1 = collist[k]; 719 v[j1] = v[j1] + d[j1] - min; 720 } 721 722 // reset row and column assignments along the alternating path. 723 do 724 { 725 i = int(pred[endofpath]); 726 colsol[endofpath] = i; 727 j1 = endofpath; 728 endofpath = rowsol[i]; 729 rowsol[i] = j1; 730 }while (i != freerow); 731 } 732 733 // calculate symmetric shape context cost 734 cv::Mat trueCostMatrix(costMatrix, cv::Rect(0,0,sizeScd1, sizeScd2)); 735 float leftcost = 0; 736 for (int nrow=0; nrow<trueCostMatrix.rows; nrow++) 737 { 738 double minval; 739 minMaxIdx(trueCostMatrix.row(nrow), &minval); 740 leftcost+=float(minval); 741 } 742 leftcost /= trueCostMatrix.rows; 743 744 float rightcost = 0; 745 for (int ncol=0; ncol<trueCostMatrix.cols; ncol++) 746 { 747 double minval; 748 minMaxIdx(trueCostMatrix.col(ncol), &minval); 749 rightcost+=float(minval); 750 } 751 rightcost /= trueCostMatrix.cols; 752 753 minMatchCost = std::max(leftcost,rightcost); 754 755 // Save in a DMatch vector 756 for (i=0;i<costMatrix.cols;i++) 757 { 758 cv::DMatch singleMatch(colsol[i],i,costMatrix.at<float>(colsol[i],i));//queryIdx,trainIdx,distance 759 outMatches.push_back(singleMatch); 760 } 761 762 // Update inliers 763 inliers1.reserve(sizeScd1); 764 for (size_t kc = 0; kc<inliers1.size(); kc++) 765 { 766 if (rowsol[kc]<sizeScd1) // if a real match 767 inliers1[kc]=1; 768 else 769 inliers1[kc]=0; 770 } 771 inliers2.reserve(sizeScd2); 772 for (size_t kc = 0; kc<inliers2.size(); kc++) 773 { 774 if (colsol[kc]<sizeScd2) // if a real match 775 inliers2[kc]=1; 776 else 777 inliers2[kc]=0; 778 } 779 } 780 781 } 782
