1 /********************************************************************* 2 * Software License Agreement (BSD License) 3 * 4 * Copyright (C) 2011 The Autonomous Systems Lab (ASL), ETH Zurich, 5 * Stefan Leutenegger, Simon Lynen and Margarita Chli. 6 * Copyright (c) 2009, Willow Garage, Inc. 7 * All rights reserved. 8 * 9 * Redistribution and use in source and binary forms, with or without 10 * modification, are permitted provided that the following conditions 11 * are met: 12 * 13 * * Redistributions of source code must retain the above copyright 14 * notice, this list of conditions and the following disclaimer. 15 * * Redistributions in binary form must reproduce the above 16 * copyright notice, this list of conditions and the following 17 * disclaimer in the documentation and/or other materials provided 18 * with the distribution. 19 * * Neither the name of the Willow Garage nor the names of its 20 * contributors may be used to endorse or promote products derived 21 * from this software without specific prior written permission. 22 * 23 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 24 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 25 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 26 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE 27 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, 28 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 29 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 30 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 31 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 32 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN 33 * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 34 * POSSIBILITY OF SUCH DAMAGE. 35 *********************************************************************/ 36 37 /* 38 BRISK - Binary Robust Invariant Scalable Keypoints 39 Reference implementation of 40 [1] Stefan Leutenegger,Margarita Chli and Roland Siegwart, BRISK: 41 Binary Robust Invariant Scalable Keypoints, in Proceedings of 42 the IEEE International Conference on Computer Vision (ICCV2011). 43 */ 44 45 #include "precomp.hpp" 46 #include <fstream> 47 #include <stdlib.h> 48 49 #include "agast_score.hpp" 50 51 namespace cv 52 { 53 54 class BRISK_Impl : public BRISK 55 { 56 public: 57 explicit BRISK_Impl(int thresh=30, int octaves=3, float patternScale=1.0f); 58 // custom setup 59 explicit BRISK_Impl(const std::vector<float> &radiusList, const std::vector<int> &numberList, 60 float dMax=5.85f, float dMin=8.2f, const std::vector<int> indexChange=std::vector<int>()); 61 62 virtual ~BRISK_Impl(); 63 64 int descriptorSize() const 65 { 66 return strings_; 67 } 68 69 int descriptorType() const 70 { 71 return CV_8U; 72 } 73 74 int defaultNorm() const 75 { 76 return NORM_HAMMING; 77 } 78 79 // call this to generate the kernel: 80 // circle of radius r (pixels), with n points; 81 // short pairings with dMax, long pairings with dMin 82 void generateKernel(const std::vector<float> &radiusList, 83 const std::vector<int> &numberList, float dMax=5.85f, float dMin=8.2f, 84 const std::vector<int> &indexChange=std::vector<int>()); 85 86 void detectAndCompute( InputArray image, InputArray mask, 87 CV_OUT std::vector<KeyPoint>& keypoints, 88 OutputArray descriptors, 89 bool useProvidedKeypoints ); 90 91 protected: 92 93 void computeKeypointsNoOrientation(InputArray image, InputArray mask, std::vector<KeyPoint>& keypoints) const; 94 void computeDescriptorsAndOrOrientation(InputArray image, InputArray mask, std::vector<KeyPoint>& keypoints, 95 OutputArray descriptors, bool doDescriptors, bool doOrientation, 96 bool useProvidedKeypoints) const; 97 98 // Feature parameters 99 CV_PROP_RW int threshold; 100 CV_PROP_RW int octaves; 101 102 // some helper structures for the Brisk pattern representation 103 struct BriskPatternPoint{ 104 float x; // x coordinate relative to center 105 float y; // x coordinate relative to center 106 float sigma; // Gaussian smoothing sigma 107 }; 108 struct BriskShortPair{ 109 unsigned int i; // index of the first pattern point 110 unsigned int j; // index of other pattern point 111 }; 112 struct BriskLongPair{ 113 unsigned int i; // index of the first pattern point 114 unsigned int j; // index of other pattern point 115 int weighted_dx; // 1024.0/dx 116 int weighted_dy; // 1024.0/dy 117 }; 118 inline int smoothedIntensity(const cv::Mat& image, 119 const cv::Mat& integral,const float key_x, 120 const float key_y, const unsigned int scale, 121 const unsigned int rot, const unsigned int point) const; 122 // pattern properties 123 BriskPatternPoint* patternPoints_; //[i][rotation][scale] 124 unsigned int points_; // total number of collocation points 125 float* scaleList_; // lists the scaling per scale index [scale] 126 unsigned int* sizeList_; // lists the total pattern size per scale index [scale] 127 static const unsigned int scales_; // scales discretization 128 static const float scalerange_; // span of sizes 40->4 Octaves - else, this needs to be adjusted... 129 static const unsigned int n_rot_; // discretization of the rotation look-up 130 131 // pairs 132 int strings_; // number of uchars the descriptor consists of 133 float dMax_; // short pair maximum distance 134 float dMin_; // long pair maximum distance 135 BriskShortPair* shortPairs_; // d<_dMax 136 BriskLongPair* longPairs_; // d>_dMin 137 unsigned int noShortPairs_; // number of shortParis 138 unsigned int noLongPairs_; // number of longParis 139 140 // general 141 static const float basicSize_; 142 }; 143 144 145 // a layer in the Brisk detector pyramid 146 class CV_EXPORTS BriskLayer 147 { 148 public: 149 // constructor arguments 150 struct CV_EXPORTS CommonParams 151 { 152 static const int HALFSAMPLE = 0; 153 static const int TWOTHIRDSAMPLE = 1; 154 }; 155 // construct a base layer 156 BriskLayer(const cv::Mat& img, float scale = 1.0f, float offset = 0.0f); 157 // derive a layer 158 BriskLayer(const BriskLayer& layer, int mode); 159 160 // Agast without non-max suppression 161 void 162 getAgastPoints(int threshold, std::vector<cv::KeyPoint>& keypoints); 163 164 // get scores - attention, this is in layer coordinates, not scale=1 coordinates! 165 inline int 166 getAgastScore(int x, int y, int threshold) const; 167 inline int 168 getAgastScore_5_8(int x, int y, int threshold) const; 169 inline int 170 getAgastScore(float xf, float yf, int threshold, float scale = 1.0f) const; 171 172 // accessors 173 inline const cv::Mat& 174 img() const 175 { 176 return img_; 177 } 178 inline const cv::Mat& 179 scores() const 180 { 181 return scores_; 182 } 183 inline float 184 scale() const 185 { 186 return scale_; 187 } 188 inline float 189 offset() const 190 { 191 return offset_; 192 } 193 194 // half sampling 195 static inline void 196 halfsample(const cv::Mat& srcimg, cv::Mat& dstimg); 197 // two third sampling 198 static inline void 199 twothirdsample(const cv::Mat& srcimg, cv::Mat& dstimg); 200 201 private: 202 // access gray values (smoothed/interpolated) 203 inline int 204 value(const cv::Mat& mat, float xf, float yf, float scale) const; 205 // the image 206 cv::Mat img_; 207 // its Agast scores 208 cv::Mat_<uchar> scores_; 209 // coordinate transformation 210 float scale_; 211 float offset_; 212 // agast 213 cv::Ptr<cv::AgastFeatureDetector> oast_9_16_; 214 int pixel_5_8_[25]; 215 int pixel_9_16_[25]; 216 }; 217 218 class CV_EXPORTS BriskScaleSpace 219 { 220 public: 221 // construct telling the octaves number: 222 BriskScaleSpace(int _octaves = 3); 223 ~BriskScaleSpace(); 224 225 // construct the image pyramids 226 void 227 constructPyramid(const cv::Mat& image); 228 229 // get Keypoints 230 void 231 getKeypoints(const int _threshold, std::vector<cv::KeyPoint>& keypoints); 232 233 protected: 234 // nonmax suppression: 235 inline bool 236 isMax2D(const int layer, const int x_layer, const int y_layer); 237 // 1D (scale axis) refinement: 238 inline float 239 refine1D(const float s_05, const float s0, const float s05, float& max) const; // around octave 240 inline float 241 refine1D_1(const float s_05, const float s0, const float s05, float& max) const; // around intra 242 inline float 243 refine1D_2(const float s_05, const float s0, const float s05, float& max) const; // around octave 0 only 244 // 2D maximum refinement: 245 inline float 246 subpixel2D(const int s_0_0, const int s_0_1, const int s_0_2, const int s_1_0, const int s_1_1, const int s_1_2, 247 const int s_2_0, const int s_2_1, const int s_2_2, float& delta_x, float& delta_y) const; 248 // 3D maximum refinement centered around (x_layer,y_layer) 249 inline float 250 refine3D(const int layer, const int x_layer, const int y_layer, float& x, float& y, float& scale, bool& ismax) const; 251 252 // interpolated score access with recalculation when needed: 253 inline int 254 getScoreAbove(const int layer, const int x_layer, const int y_layer) const; 255 inline int 256 getScoreBelow(const int layer, const int x_layer, const int y_layer) const; 257 258 // return the maximum of score patches above or below 259 inline float 260 getScoreMaxAbove(const int layer, const int x_layer, const int y_layer, const int threshold, bool& ismax, 261 float& dx, float& dy) const; 262 inline float 263 getScoreMaxBelow(const int layer, const int x_layer, const int y_layer, const int threshold, bool& ismax, 264 float& dx, float& dy) const; 265 266 // the image pyramids: 267 int layers_; 268 std::vector<BriskLayer> pyramid_; 269 270 // some constant parameters: 271 static const float safetyFactor_; 272 static const float basicSize_; 273 }; 274 275 const float BRISK_Impl::basicSize_ = 12.0f; 276 const unsigned int BRISK_Impl::scales_ = 64; 277 const float BRISK_Impl::scalerange_ = 30.f; // 40->4 Octaves - else, this needs to be adjusted... 278 const unsigned int BRISK_Impl::n_rot_ = 1024; // discretization of the rotation look-up 279 280 const float BriskScaleSpace::safetyFactor_ = 1.0f; 281 const float BriskScaleSpace::basicSize_ = 12.0f; 282 283 // constructors 284 BRISK_Impl::BRISK_Impl(int thresh, int octaves_in, float patternScale) 285 { 286 threshold = thresh; 287 octaves = octaves_in; 288 289 std::vector<float> rList; 290 std::vector<int> nList; 291 292 // this is the standard pattern found to be suitable also 293 rList.resize(5); 294 nList.resize(5); 295 const double f = 0.85 * patternScale; 296 297 rList[0] = (float)(f * 0.); 298 rList[1] = (float)(f * 2.9); 299 rList[2] = (float)(f * 4.9); 300 rList[3] = (float)(f * 7.4); 301 rList[4] = (float)(f * 10.8); 302 303 nList[0] = 1; 304 nList[1] = 10; 305 nList[2] = 14; 306 nList[3] = 15; 307 nList[4] = 20; 308 309 generateKernel(rList, nList, (float)(5.85 * patternScale), (float)(8.2 * patternScale)); 310 } 311 312 BRISK_Impl::BRISK_Impl(const std::vector<float> &radiusList, 313 const std::vector<int> &numberList, 314 float dMax, float dMin, 315 const std::vector<int> indexChange) 316 { 317 generateKernel(radiusList, numberList, dMax, dMin, indexChange); 318 threshold = 20; 319 octaves = 3; 320 } 321 322 void 323 BRISK_Impl::generateKernel(const std::vector<float> &radiusList, 324 const std::vector<int> &numberList, 325 float dMax, float dMin, 326 const std::vector<int>& _indexChange) 327 { 328 std::vector<int> indexChange = _indexChange; 329 dMax_ = dMax; 330 dMin_ = dMin; 331 332 // get the total number of points 333 const int rings = (int)radiusList.size(); 334 CV_Assert(radiusList.size() != 0 && radiusList.size() == numberList.size()); 335 points_ = 0; // remember the total number of points 336 for (int ring = 0; ring < rings; ring++) 337 { 338 points_ += numberList[ring]; 339 } 340 // set up the patterns 341 patternPoints_ = new BriskPatternPoint[points_ * scales_ * n_rot_]; 342 BriskPatternPoint* patternIterator = patternPoints_; 343 344 // define the scale discretization: 345 static const float lb_scale = (float)(std::log(scalerange_) / std::log(2.0)); 346 static const float lb_scale_step = lb_scale / (scales_); 347 348 scaleList_ = new float[scales_]; 349 sizeList_ = new unsigned int[scales_]; 350 351 const float sigma_scale = 1.3f; 352 353 for (unsigned int scale = 0; scale < scales_; ++scale) 354 { 355 scaleList_[scale] = (float)std::pow((double) 2.0, (double) (scale * lb_scale_step)); 356 sizeList_[scale] = 0; 357 358 // generate the pattern points look-up 359 double alpha, theta; 360 for (size_t rot = 0; rot < n_rot_; ++rot) 361 { 362 theta = double(rot) * 2 * CV_PI / double(n_rot_); // this is the rotation of the feature 363 for (int ring = 0; ring < rings; ++ring) 364 { 365 for (int num = 0; num < numberList[ring]; ++num) 366 { 367 // the actual coordinates on the circle 368 alpha = (double(num)) * 2 * CV_PI / double(numberList[ring]); 369 patternIterator->x = (float)(scaleList_[scale] * radiusList[ring] * cos(alpha + theta)); // feature rotation plus angle of the point 370 patternIterator->y = (float)(scaleList_[scale] * radiusList[ring] * sin(alpha + theta)); 371 // and the gaussian kernel sigma 372 if (ring == 0) 373 { 374 patternIterator->sigma = sigma_scale * scaleList_[scale] * 0.5f; 375 } 376 else 377 { 378 patternIterator->sigma = (float)(sigma_scale * scaleList_[scale] * (double(radiusList[ring])) 379 * sin(CV_PI / numberList[ring])); 380 } 381 // adapt the sizeList if necessary 382 const unsigned int size = cvCeil(((scaleList_[scale] * radiusList[ring]) + patternIterator->sigma)) + 1; 383 if (sizeList_[scale] < size) 384 { 385 sizeList_[scale] = size; 386 } 387 388 // increment the iterator 389 ++patternIterator; 390 } 391 } 392 } 393 } 394 395 // now also generate pairings 396 shortPairs_ = new BriskShortPair[points_ * (points_ - 1) / 2]; 397 longPairs_ = new BriskLongPair[points_ * (points_ - 1) / 2]; 398 noShortPairs_ = 0; 399 noLongPairs_ = 0; 400 401 // fill indexChange with 0..n if empty 402 unsigned int indSize = (unsigned int)indexChange.size(); 403 if (indSize == 0) 404 { 405 indexChange.resize(points_ * (points_ - 1) / 2); 406 indSize = (unsigned int)indexChange.size(); 407 408 for (unsigned int i = 0; i < indSize; i++) 409 indexChange[i] = i; 410 } 411 const float dMin_sq = dMin_ * dMin_; 412 const float dMax_sq = dMax_ * dMax_; 413 for (unsigned int i = 1; i < points_; i++) 414 { 415 for (unsigned int j = 0; j < i; j++) 416 { //(find all the pairs) 417 // point pair distance: 418 const float dx = patternPoints_[j].x - patternPoints_[i].x; 419 const float dy = patternPoints_[j].y - patternPoints_[i].y; 420 const float norm_sq = (dx * dx + dy * dy); 421 if (norm_sq > dMin_sq) 422 { 423 // save to long pairs 424 BriskLongPair& longPair = longPairs_[noLongPairs_]; 425 longPair.weighted_dx = int((dx / (norm_sq)) * 2048.0 + 0.5); 426 longPair.weighted_dy = int((dy / (norm_sq)) * 2048.0 + 0.5); 427 longPair.i = i; 428 longPair.j = j; 429 ++noLongPairs_; 430 } 431 else if (norm_sq < dMax_sq) 432 { 433 // save to short pairs 434 CV_Assert(noShortPairs_ < indSize); 435 // make sure the user passes something sensible 436 BriskShortPair& shortPair = shortPairs_[indexChange[noShortPairs_]]; 437 shortPair.j = j; 438 shortPair.i = i; 439 ++noShortPairs_; 440 } 441 } 442 } 443 444 // no bits: 445 strings_ = (int) ceil((float(noShortPairs_)) / 128.0) * 4 * 4; 446 } 447 448 // simple alternative: 449 inline int 450 BRISK_Impl::smoothedIntensity(const cv::Mat& image, const cv::Mat& integral, const float key_x, 451 const float key_y, const unsigned int scale, const unsigned int rot, 452 const unsigned int point) const 453 { 454 455 // get the float position 456 const BriskPatternPoint& briskPoint = patternPoints_[scale * n_rot_ * points_ + rot * points_ + point]; 457 const float xf = briskPoint.x + key_x; 458 const float yf = briskPoint.y + key_y; 459 const int x = int(xf); 460 const int y = int(yf); 461 const int& imagecols = image.cols; 462 463 // get the sigma: 464 const float sigma_half = briskPoint.sigma; 465 const float area = 4.0f * sigma_half * sigma_half; 466 467 // calculate output: 468 int ret_val; 469 if (sigma_half < 0.5) 470 { 471 //interpolation multipliers: 472 const int r_x = (int)((xf - x) * 1024); 473 const int r_y = (int)((yf - y) * 1024); 474 const int r_x_1 = (1024 - r_x); 475 const int r_y_1 = (1024 - r_y); 476 const uchar* ptr = &image.at<uchar>(y, x); 477 size_t step = image.step; 478 // just interpolate: 479 ret_val = r_x_1 * r_y_1 * ptr[0] + r_x * r_y_1 * ptr[1] + 480 r_x * r_y * ptr[step] + r_x_1 * r_y * ptr[step+1]; 481 return (ret_val + 512) / 1024; 482 } 483 484 // this is the standard case (simple, not speed optimized yet): 485 486 // scaling: 487 const int scaling = (int)(4194304.0 / area); 488 const int scaling2 = int(float(scaling) * area / 1024.0); 489 490 // the integral image is larger: 491 const int integralcols = imagecols + 1; 492 493 // calculate borders 494 const float x_1 = xf - sigma_half; 495 const float x1 = xf + sigma_half; 496 const float y_1 = yf - sigma_half; 497 const float y1 = yf + sigma_half; 498 499 const int x_left = int(x_1 + 0.5); 500 const int y_top = int(y_1 + 0.5); 501 const int x_right = int(x1 + 0.5); 502 const int y_bottom = int(y1 + 0.5); 503 504 // overlap area - multiplication factors: 505 const float r_x_1 = float(x_left) - x_1 + 0.5f; 506 const float r_y_1 = float(y_top) - y_1 + 0.5f; 507 const float r_x1 = x1 - float(x_right) + 0.5f; 508 const float r_y1 = y1 - float(y_bottom) + 0.5f; 509 const int dx = x_right - x_left - 1; 510 const int dy = y_bottom - y_top - 1; 511 const int A = (int)((r_x_1 * r_y_1) * scaling); 512 const int B = (int)((r_x1 * r_y_1) * scaling); 513 const int C = (int)((r_x1 * r_y1) * scaling); 514 const int D = (int)((r_x_1 * r_y1) * scaling); 515 const int r_x_1_i = (int)(r_x_1 * scaling); 516 const int r_y_1_i = (int)(r_y_1 * scaling); 517 const int r_x1_i = (int)(r_x1 * scaling); 518 const int r_y1_i = (int)(r_y1 * scaling); 519 520 if (dx + dy > 2) 521 { 522 // now the calculation: 523 const uchar* ptr = image.ptr() + x_left + imagecols * y_top; 524 // first the corners: 525 ret_val = A * int(*ptr); 526 ptr += dx + 1; 527 ret_val += B * int(*ptr); 528 ptr += dy * imagecols + 1; 529 ret_val += C * int(*ptr); 530 ptr -= dx + 1; 531 ret_val += D * int(*ptr); 532 533 // next the edges: 534 const int* ptr_integral = integral.ptr<int>() + x_left + integralcols * y_top + 1; 535 // find a simple path through the different surface corners 536 const int tmp1 = (*ptr_integral); 537 ptr_integral += dx; 538 const int tmp2 = (*ptr_integral); 539 ptr_integral += integralcols; 540 const int tmp3 = (*ptr_integral); 541 ptr_integral++; 542 const int tmp4 = (*ptr_integral); 543 ptr_integral += dy * integralcols; 544 const int tmp5 = (*ptr_integral); 545 ptr_integral--; 546 const int tmp6 = (*ptr_integral); 547 ptr_integral += integralcols; 548 const int tmp7 = (*ptr_integral); 549 ptr_integral -= dx; 550 const int tmp8 = (*ptr_integral); 551 ptr_integral -= integralcols; 552 const int tmp9 = (*ptr_integral); 553 ptr_integral--; 554 const int tmp10 = (*ptr_integral); 555 ptr_integral -= dy * integralcols; 556 const int tmp11 = (*ptr_integral); 557 ptr_integral++; 558 const int tmp12 = (*ptr_integral); 559 560 // assign the weighted surface integrals: 561 const int upper = (tmp3 - tmp2 + tmp1 - tmp12) * r_y_1_i; 562 const int middle = (tmp6 - tmp3 + tmp12 - tmp9) * scaling; 563 const int left = (tmp9 - tmp12 + tmp11 - tmp10) * r_x_1_i; 564 const int right = (tmp5 - tmp4 + tmp3 - tmp6) * r_x1_i; 565 const int bottom = (tmp7 - tmp6 + tmp9 - tmp8) * r_y1_i; 566 567 return (ret_val + upper + middle + left + right + bottom + scaling2 / 2) / scaling2; 568 } 569 570 // now the calculation: 571 const uchar* ptr = image.ptr() + x_left + imagecols * y_top; 572 // first row: 573 ret_val = A * int(*ptr); 574 ptr++; 575 const uchar* end1 = ptr + dx; 576 for (; ptr < end1; ptr++) 577 { 578 ret_val += r_y_1_i * int(*ptr); 579 } 580 ret_val += B * int(*ptr); 581 // middle ones: 582 ptr += imagecols - dx - 1; 583 const uchar* end_j = ptr + dy * imagecols; 584 for (; ptr < end_j; ptr += imagecols - dx - 1) 585 { 586 ret_val += r_x_1_i * int(*ptr); 587 ptr++; 588 const uchar* end2 = ptr + dx; 589 for (; ptr < end2; ptr++) 590 { 591 ret_val += int(*ptr) * scaling; 592 } 593 ret_val += r_x1_i * int(*ptr); 594 } 595 // last row: 596 ret_val += D * int(*ptr); 597 ptr++; 598 const uchar* end3 = ptr + dx; 599 for (; ptr < end3; ptr++) 600 { 601 ret_val += r_y1_i * int(*ptr); 602 } 603 ret_val += C * int(*ptr); 604 605 return (ret_val + scaling2 / 2) / scaling2; 606 } 607 608 inline bool 609 RoiPredicate(const float minX, const float minY, const float maxX, const float maxY, const KeyPoint& keyPt) 610 { 611 const Point2f& pt = keyPt.pt; 612 return (pt.x < minX) || (pt.x >= maxX) || (pt.y < minY) || (pt.y >= maxY); 613 } 614 615 // computes the descriptor 616 void 617 BRISK_Impl::detectAndCompute( InputArray _image, InputArray _mask, std::vector<KeyPoint>& keypoints, 618 OutputArray _descriptors, bool useProvidedKeypoints) 619 { 620 bool doOrientation=true; 621 622 // If the user specified cv::noArray(), this will yield false. Otherwise it will return true. 623 bool doDescriptors = _descriptors.needed(); 624 625 computeDescriptorsAndOrOrientation(_image, _mask, keypoints, _descriptors, doDescriptors, doOrientation, 626 useProvidedKeypoints); 627 } 628 629 void 630 BRISK_Impl::computeDescriptorsAndOrOrientation(InputArray _image, InputArray _mask, std::vector<KeyPoint>& keypoints, 631 OutputArray _descriptors, bool doDescriptors, bool doOrientation, 632 bool useProvidedKeypoints) const 633 { 634 Mat image = _image.getMat(), mask = _mask.getMat(); 635 if( image.type() != CV_8UC1 ) 636 cvtColor(image, image, COLOR_BGR2GRAY); 637 638 if (!useProvidedKeypoints) 639 { 640 doOrientation = true; 641 computeKeypointsNoOrientation(_image, _mask, keypoints); 642 } 643 644 //Remove keypoints very close to the border 645 size_t ksize = keypoints.size(); 646 std::vector<int> kscales; // remember the scale per keypoint 647 kscales.resize(ksize); 648 static const float log2 = 0.693147180559945f; 649 static const float lb_scalerange = (float)(std::log(scalerange_) / (log2)); 650 std::vector<cv::KeyPoint>::iterator beginning = keypoints.begin(); 651 std::vector<int>::iterator beginningkscales = kscales.begin(); 652 static const float basicSize06 = basicSize_ * 0.6f; 653 for (size_t k = 0; k < ksize; k++) 654 { 655 unsigned int scale; 656 scale = std::max((int) (scales_ / lb_scalerange * (std::log(keypoints[k].size / (basicSize06)) / log2) + 0.5), 0); 657 // saturate 658 if (scale >= scales_) 659 scale = scales_ - 1; 660 kscales[k] = scale; 661 const int border = sizeList_[scale]; 662 const int border_x = image.cols - border; 663 const int border_y = image.rows - border; 664 if (RoiPredicate((float)border, (float)border, (float)border_x, (float)border_y, keypoints[k])) 665 { 666 keypoints.erase(beginning + k); 667 kscales.erase(beginningkscales + k); 668 if (k == 0) 669 { 670 beginning = keypoints.begin(); 671 beginningkscales = kscales.begin(); 672 } 673 ksize--; 674 k--; 675 } 676 } 677 678 // first, calculate the integral image over the whole image: 679 // current integral image 680 cv::Mat _integral; // the integral image 681 cv::integral(image, _integral); 682 683 int* _values = new int[points_]; // for temporary use 684 685 // resize the descriptors: 686 cv::Mat descriptors; 687 if (doDescriptors) 688 { 689 _descriptors.create((int)ksize, strings_, CV_8U); 690 descriptors = _descriptors.getMat(); 691 descriptors.setTo(0); 692 } 693 694 // now do the extraction for all keypoints: 695 696 // temporary variables containing gray values at sample points: 697 int t1; 698 int t2; 699 700 // the feature orientation 701 const uchar* ptr = descriptors.ptr(); 702 for (size_t k = 0; k < ksize; k++) 703 { 704 cv::KeyPoint& kp = keypoints[k]; 705 const int& scale = kscales[k]; 706 int* pvalues = _values; 707 const float& x = kp.pt.x; 708 const float& y = kp.pt.y; 709 710 if (doOrientation) 711 { 712 // get the gray values in the unrotated pattern 713 for (unsigned int i = 0; i < points_; i++) 714 { 715 *(pvalues++) = smoothedIntensity(image, _integral, x, y, scale, 0, i); 716 } 717 718 int direction0 = 0; 719 int direction1 = 0; 720 // now iterate through the long pairings 721 const BriskLongPair* max = longPairs_ + noLongPairs_; 722 for (BriskLongPair* iter = longPairs_; iter < max; ++iter) 723 { 724 t1 = *(_values + iter->i); 725 t2 = *(_values + iter->j); 726 const int delta_t = (t1 - t2); 727 // update the direction: 728 const int tmp0 = delta_t * (iter->weighted_dx) / 1024; 729 const int tmp1 = delta_t * (iter->weighted_dy) / 1024; 730 direction0 += tmp0; 731 direction1 += tmp1; 732 } 733 kp.angle = (float)(atan2((float) direction1, (float) direction0) / CV_PI * 180.0); 734 735 if (!doDescriptors) 736 { 737 if (kp.angle < 0) 738 kp.angle += 360.f; 739 } 740 } 741 742 if (!doDescriptors) 743 continue; 744 745 int theta; 746 if (kp.angle==-1) 747 { 748 // don't compute the gradient direction, just assign a rotation of 0 749 theta = 0; 750 } 751 else 752 { 753 theta = (int) (n_rot_ * (kp.angle / (360.0)) + 0.5); 754 if (theta < 0) 755 theta += n_rot_; 756 if (theta >= int(n_rot_)) 757 theta -= n_rot_; 758 } 759 760 if (kp.angle < 0) 761 kp.angle += 360.f; 762 763 // now also extract the stuff for the actual direction: 764 // let us compute the smoothed values 765 int shifter = 0; 766 767 //unsigned int mean=0; 768 pvalues = _values; 769 // get the gray values in the rotated pattern 770 for (unsigned int i = 0; i < points_; i++) 771 { 772 *(pvalues++) = smoothedIntensity(image, _integral, x, y, scale, theta, i); 773 } 774 775 // now iterate through all the pairings 776 unsigned int* ptr2 = (unsigned int*) ptr; 777 const BriskShortPair* max = shortPairs_ + noShortPairs_; 778 for (BriskShortPair* iter = shortPairs_; iter < max; ++iter) 779 { 780 t1 = *(_values + iter->i); 781 t2 = *(_values + iter->j); 782 if (t1 > t2) 783 { 784 *ptr2 |= ((1) << shifter); 785 786 } // else already initialized with zero 787 // take care of the iterators: 788 ++shifter; 789 if (shifter == 32) 790 { 791 shifter = 0; 792 ++ptr2; 793 } 794 } 795 796 ptr += strings_; 797 } 798 799 // clean-up 800 delete[] _values; 801 } 802 803 804 BRISK_Impl::~BRISK_Impl() 805 { 806 delete[] patternPoints_; 807 delete[] shortPairs_; 808 delete[] longPairs_; 809 delete[] scaleList_; 810 delete[] sizeList_; 811 } 812 813 void 814 BRISK_Impl::computeKeypointsNoOrientation(InputArray _image, InputArray _mask, std::vector<KeyPoint>& keypoints) const 815 { 816 Mat image = _image.getMat(), mask = _mask.getMat(); 817 if( image.type() != CV_8UC1 ) 818 cvtColor(_image, image, COLOR_BGR2GRAY); 819 820 BriskScaleSpace briskScaleSpace(octaves); 821 briskScaleSpace.constructPyramid(image); 822 briskScaleSpace.getKeypoints(threshold, keypoints); 823 824 // remove invalid points 825 KeyPointsFilter::runByPixelsMask(keypoints, mask); 826 } 827 828 // construct telling the octaves number: 829 BriskScaleSpace::BriskScaleSpace(int _octaves) 830 { 831 if (_octaves == 0) 832 layers_ = 1; 833 else 834 layers_ = 2 * _octaves; 835 } 836 BriskScaleSpace::~BriskScaleSpace() 837 { 838 839 } 840 // construct the image pyramids 841 void 842 BriskScaleSpace::constructPyramid(const cv::Mat& image) 843 { 844 845 // set correct size: 846 pyramid_.clear(); 847 848 // fill the pyramid: 849 pyramid_.push_back(BriskLayer(image.clone())); 850 if (layers_ > 1) 851 { 852 pyramid_.push_back(BriskLayer(pyramid_.back(), BriskLayer::CommonParams::TWOTHIRDSAMPLE)); 853 } 854 const int octaves2 = layers_; 855 856 for (uchar i = 2; i < octaves2; i += 2) 857 { 858 pyramid_.push_back(BriskLayer(pyramid_[i - 2], BriskLayer::CommonParams::HALFSAMPLE)); 859 pyramid_.push_back(BriskLayer(pyramid_[i - 1], BriskLayer::CommonParams::HALFSAMPLE)); 860 } 861 } 862 863 void 864 BriskScaleSpace::getKeypoints(const int threshold_, std::vector<cv::KeyPoint>& keypoints) 865 { 866 // make sure keypoints is empty 867 keypoints.resize(0); 868 keypoints.reserve(2000); 869 870 // assign thresholds 871 int safeThreshold_ = (int)(threshold_ * safetyFactor_); 872 std::vector<std::vector<cv::KeyPoint> > agastPoints; 873 agastPoints.resize(layers_); 874 875 // go through the octaves and intra layers and calculate agast corner scores: 876 for (int i = 0; i < layers_; i++) 877 { 878 // call OAST16_9 without nms 879 BriskLayer& l = pyramid_[i]; 880 l.getAgastPoints(safeThreshold_, agastPoints[i]); 881 } 882 883 if (layers_ == 1) 884 { 885 // just do a simple 2d subpixel refinement... 886 const size_t num = agastPoints[0].size(); 887 for (size_t n = 0; n < num; n++) 888 { 889 const cv::Point2f& point = agastPoints.at(0)[n].pt; 890 // first check if it is a maximum: 891 if (!isMax2D(0, (int)point.x, (int)point.y)) 892 continue; 893 894 // let's do the subpixel and float scale refinement: 895 BriskLayer& l = pyramid_[0]; 896 int s_0_0 = l.getAgastScore(point.x - 1, point.y - 1, 1); 897 int s_1_0 = l.getAgastScore(point.x, point.y - 1, 1); 898 int s_2_0 = l.getAgastScore(point.x + 1, point.y - 1, 1); 899 int s_2_1 = l.getAgastScore(point.x + 1, point.y, 1); 900 int s_1_1 = l.getAgastScore(point.x, point.y, 1); 901 int s_0_1 = l.getAgastScore(point.x - 1, point.y, 1); 902 int s_0_2 = l.getAgastScore(point.x - 1, point.y + 1, 1); 903 int s_1_2 = l.getAgastScore(point.x, point.y + 1, 1); 904 int s_2_2 = l.getAgastScore(point.x + 1, point.y + 1, 1); 905 float delta_x, delta_y; 906 float max = subpixel2D(s_0_0, s_0_1, s_0_2, s_1_0, s_1_1, s_1_2, s_2_0, s_2_1, s_2_2, delta_x, delta_y); 907 908 // store: 909 keypoints.push_back(cv::KeyPoint(float(point.x) + delta_x, float(point.y) + delta_y, basicSize_, -1, max, 0)); 910 911 } 912 913 return; 914 } 915 916 float x, y, scale, score; 917 for (int i = 0; i < layers_; i++) 918 { 919 BriskLayer& l = pyramid_[i]; 920 const size_t num = agastPoints[i].size(); 921 if (i == layers_ - 1) 922 { 923 for (size_t n = 0; n < num; n++) 924 { 925 const cv::Point2f& point = agastPoints.at(i)[n].pt; 926 // consider only 2D maxima... 927 if (!isMax2D(i, (int)point.x, (int)point.y)) 928 continue; 929 930 bool ismax; 931 float dx, dy; 932 getScoreMaxBelow(i, (int)point.x, (int)point.y, l.getAgastScore(point.x, point.y, safeThreshold_), ismax, dx, dy); 933 if (!ismax) 934 continue; 935 936 // get the patch on this layer: 937 int s_0_0 = l.getAgastScore(point.x - 1, point.y - 1, 1); 938 int s_1_0 = l.getAgastScore(point.x, point.y - 1, 1); 939 int s_2_0 = l.getAgastScore(point.x + 1, point.y - 1, 1); 940 int s_2_1 = l.getAgastScore(point.x + 1, point.y, 1); 941 int s_1_1 = l.getAgastScore(point.x, point.y, 1); 942 int s_0_1 = l.getAgastScore(point.x - 1, point.y, 1); 943 int s_0_2 = l.getAgastScore(point.x - 1, point.y + 1, 1); 944 int s_1_2 = l.getAgastScore(point.x, point.y + 1, 1); 945 int s_2_2 = l.getAgastScore(point.x + 1, point.y + 1, 1); 946 float delta_x, delta_y; 947 float max = subpixel2D(s_0_0, s_0_1, s_0_2, s_1_0, s_1_1, s_1_2, s_2_0, s_2_1, s_2_2, delta_x, delta_y); 948 949 // store: 950 keypoints.push_back( 951 cv::KeyPoint((float(point.x) + delta_x) * l.scale() + l.offset(), 952 (float(point.y) + delta_y) * l.scale() + l.offset(), basicSize_ * l.scale(), -1, max, i)); 953 } 954 } 955 else 956 { 957 // not the last layer: 958 for (size_t n = 0; n < num; n++) 959 { 960 const cv::Point2f& point = agastPoints.at(i)[n].pt; 961 962 // first check if it is a maximum: 963 if (!isMax2D(i, (int)point.x, (int)point.y)) 964 continue; 965 966 // let's do the subpixel and float scale refinement: 967 bool ismax=false; 968 score = refine3D(i, (int)point.x, (int)point.y, x, y, scale, ismax); 969 if (!ismax) 970 { 971 continue; 972 } 973 974 // finally store the detected keypoint: 975 if (score > float(threshold_)) 976 { 977 keypoints.push_back(cv::KeyPoint(x, y, basicSize_ * scale, -1, score, i)); 978 } 979 } 980 } 981 } 982 } 983 984 // interpolated score access with recalculation when needed: 985 inline int 986 BriskScaleSpace::getScoreAbove(const int layer, const int x_layer, const int y_layer) const 987 { 988 CV_Assert(layer < layers_-1); 989 const BriskLayer& l = pyramid_[layer + 1]; 990 if (layer % 2 == 0) 991 { // octave 992 const int sixths_x = 4 * x_layer - 1; 993 const int x_above = sixths_x / 6; 994 const int sixths_y = 4 * y_layer - 1; 995 const int y_above = sixths_y / 6; 996 const int r_x = (sixths_x % 6); 997 const int r_x_1 = 6 - r_x; 998 const int r_y = (sixths_y % 6); 999 const int r_y_1 = 6 - r_y; 1000 uchar score = 0xFF 1001 & ((r_x_1 * r_y_1 * l.getAgastScore(x_above, y_above, 1) + r_x * r_y_1 1002 * l.getAgastScore(x_above + 1, y_above, 1) 1003 + r_x_1 * r_y * l.getAgastScore(x_above, y_above + 1, 1) 1004 + r_x * r_y * l.getAgastScore(x_above + 1, y_above + 1, 1) + 18) 1005 / 36); 1006 1007 return score; 1008 } 1009 else 1010 { // intra 1011 const int eighths_x = 6 * x_layer - 1; 1012 const int x_above = eighths_x / 8; 1013 const int eighths_y = 6 * y_layer - 1; 1014 const int y_above = eighths_y / 8; 1015 const int r_x = (eighths_x % 8); 1016 const int r_x_1 = 8 - r_x; 1017 const int r_y = (eighths_y % 8); 1018 const int r_y_1 = 8 - r_y; 1019 uchar score = 0xFF 1020 & ((r_x_1 * r_y_1 * l.getAgastScore(x_above, y_above, 1) + r_x * r_y_1 1021 * l.getAgastScore(x_above + 1, y_above, 1) 1022 + r_x_1 * r_y * l.getAgastScore(x_above, y_above + 1, 1) 1023 + r_x * r_y * l.getAgastScore(x_above + 1, y_above + 1, 1) + 32) 1024 / 64); 1025 return score; 1026 } 1027 } 1028 inline int 1029 BriskScaleSpace::getScoreBelow(const int layer, const int x_layer, const int y_layer) const 1030 { 1031 CV_Assert(layer); 1032 const BriskLayer& l = pyramid_[layer - 1]; 1033 int sixth_x; 1034 int quarter_x; 1035 float xf; 1036 int sixth_y; 1037 int quarter_y; 1038 float yf; 1039 1040 // scaling: 1041 float offs; 1042 float area; 1043 int scaling; 1044 int scaling2; 1045 1046 if (layer % 2 == 0) 1047 { // octave 1048 sixth_x = 8 * x_layer + 1; 1049 xf = float(sixth_x) / 6.0f; 1050 sixth_y = 8 * y_layer + 1; 1051 yf = float(sixth_y) / 6.0f; 1052 1053 // scaling: 1054 offs = 2.0f / 3.0f; 1055 area = 4.0f * offs * offs; 1056 scaling = (int)(4194304.0 / area); 1057 scaling2 = (int)(float(scaling) * area); 1058 } 1059 else 1060 { 1061 quarter_x = 6 * x_layer + 1; 1062 xf = float(quarter_x) / 4.0f; 1063 quarter_y = 6 * y_layer + 1; 1064 yf = float(quarter_y) / 4.0f; 1065 1066 // scaling: 1067 offs = 3.0f / 4.0f; 1068 area = 4.0f * offs * offs; 1069 scaling = (int)(4194304.0 / area); 1070 scaling2 = (int)(float(scaling) * area); 1071 } 1072 1073 // calculate borders 1074 const float x_1 = xf - offs; 1075 const float x1 = xf + offs; 1076 const float y_1 = yf - offs; 1077 const float y1 = yf + offs; 1078 1079 const int x_left = int(x_1 + 0.5); 1080 const int y_top = int(y_1 + 0.5); 1081 const int x_right = int(x1 + 0.5); 1082 const int y_bottom = int(y1 + 0.5); 1083 1084 // overlap area - multiplication factors: 1085 const float r_x_1 = float(x_left) - x_1 + 0.5f; 1086 const float r_y_1 = float(y_top) - y_1 + 0.5f; 1087 const float r_x1 = x1 - float(x_right) + 0.5f; 1088 const float r_y1 = y1 - float(y_bottom) + 0.5f; 1089 const int dx = x_right - x_left - 1; 1090 const int dy = y_bottom - y_top - 1; 1091 const int A = (int)((r_x_1 * r_y_1) * scaling); 1092 const int B = (int)((r_x1 * r_y_1) * scaling); 1093 const int C = (int)((r_x1 * r_y1) * scaling); 1094 const int D = (int)((r_x_1 * r_y1) * scaling); 1095 const int r_x_1_i = (int)(r_x_1 * scaling); 1096 const int r_y_1_i = (int)(r_y_1 * scaling); 1097 const int r_x1_i = (int)(r_x1 * scaling); 1098 const int r_y1_i = (int)(r_y1 * scaling); 1099 1100 // first row: 1101 int ret_val = A * int(l.getAgastScore(x_left, y_top, 1)); 1102 for (int X = 1; X <= dx; X++) 1103 { 1104 ret_val += r_y_1_i * int(l.getAgastScore(x_left + X, y_top, 1)); 1105 } 1106 ret_val += B * int(l.getAgastScore(x_left + dx + 1, y_top, 1)); 1107 // middle ones: 1108 for (int Y = 1; Y <= dy; Y++) 1109 { 1110 ret_val += r_x_1_i * int(l.getAgastScore(x_left, y_top + Y, 1)); 1111 1112 for (int X = 1; X <= dx; X++) 1113 { 1114 ret_val += int(l.getAgastScore(x_left + X, y_top + Y, 1)) * scaling; 1115 } 1116 ret_val += r_x1_i * int(l.getAgastScore(x_left + dx + 1, y_top + Y, 1)); 1117 } 1118 // last row: 1119 ret_val += D * int(l.getAgastScore(x_left, y_top + dy + 1, 1)); 1120 for (int X = 1; X <= dx; X++) 1121 { 1122 ret_val += r_y1_i * int(l.getAgastScore(x_left + X, y_top + dy + 1, 1)); 1123 } 1124 ret_val += C * int(l.getAgastScore(x_left + dx + 1, y_top + dy + 1, 1)); 1125 1126 return ((ret_val + scaling2 / 2) / scaling2); 1127 } 1128 1129 inline bool 1130 BriskScaleSpace::isMax2D(const int layer, const int x_layer, const int y_layer) 1131 { 1132 const cv::Mat& scores = pyramid_[layer].scores(); 1133 const int scorescols = scores.cols; 1134 const uchar* data = scores.ptr() + y_layer * scorescols + x_layer; 1135 // decision tree: 1136 const uchar center = (*data); 1137 data--; 1138 const uchar s_10 = *data; 1139 if (center < s_10) 1140 return false; 1141 data += 2; 1142 const uchar s10 = *data; 1143 if (center < s10) 1144 return false; 1145 data -= (scorescols + 1); 1146 const uchar s0_1 = *data; 1147 if (center < s0_1) 1148 return false; 1149 data += 2 * scorescols; 1150 const uchar s01 = *data; 1151 if (center < s01) 1152 return false; 1153 data--; 1154 const uchar s_11 = *data; 1155 if (center < s_11) 1156 return false; 1157 data += 2; 1158 const uchar s11 = *data; 1159 if (center < s11) 1160 return false; 1161 data -= 2 * scorescols; 1162 const uchar s1_1 = *data; 1163 if (center < s1_1) 1164 return false; 1165 data -= 2; 1166 const uchar s_1_1 = *data; 1167 if (center < s_1_1) 1168 return false; 1169 1170 // reject neighbor maxima 1171 std::vector<int> delta; 1172 // put together a list of 2d-offsets to where the maximum is also reached 1173 if (center == s_1_1) 1174 { 1175 delta.push_back(-1); 1176 delta.push_back(-1); 1177 } 1178 if (center == s0_1) 1179 { 1180 delta.push_back(0); 1181 delta.push_back(-1); 1182 } 1183 if (center == s1_1) 1184 { 1185 delta.push_back(1); 1186 delta.push_back(-1); 1187 } 1188 if (center == s_10) 1189 { 1190 delta.push_back(-1); 1191 delta.push_back(0); 1192 } 1193 if (center == s10) 1194 { 1195 delta.push_back(1); 1196 delta.push_back(0); 1197 } 1198 if (center == s_11) 1199 { 1200 delta.push_back(-1); 1201 delta.push_back(1); 1202 } 1203 if (center == s01) 1204 { 1205 delta.push_back(0); 1206 delta.push_back(1); 1207 } 1208 if (center == s11) 1209 { 1210 delta.push_back(1); 1211 delta.push_back(1); 1212 } 1213 const unsigned int deltasize = (unsigned int)delta.size(); 1214 if (deltasize != 0) 1215 { 1216 // in this case, we have to analyze the situation more carefully: 1217 // the values are gaussian blurred and then we really decide 1218 data = scores.ptr() + y_layer * scorescols + x_layer; 1219 int smoothedcenter = 4 * center + 2 * (s_10 + s10 + s0_1 + s01) + s_1_1 + s1_1 + s_11 + s11; 1220 for (unsigned int i = 0; i < deltasize; i += 2) 1221 { 1222 data = scores.ptr() + (y_layer - 1 + delta[i + 1]) * scorescols + x_layer + delta[i] - 1; 1223 int othercenter = *data; 1224 data++; 1225 othercenter += 2 * (*data); 1226 data++; 1227 othercenter += *data; 1228 data += scorescols; 1229 othercenter += 2 * (*data); 1230 data--; 1231 othercenter += 4 * (*data); 1232 data--; 1233 othercenter += 2 * (*data); 1234 data += scorescols; 1235 othercenter += *data; 1236 data++; 1237 othercenter += 2 * (*data); 1238 data++; 1239 othercenter += *data; 1240 if (othercenter > smoothedcenter) 1241 return false; 1242 } 1243 } 1244 return true; 1245 } 1246 1247 // 3D maximum refinement centered around (x_layer,y_layer) 1248 inline float 1249 BriskScaleSpace::refine3D(const int layer, const int x_layer, const int y_layer, float& x, float& y, float& scale, 1250 bool& ismax) const 1251 { 1252 ismax = true; 1253 const BriskLayer& thisLayer = pyramid_[layer]; 1254 const int center = thisLayer.getAgastScore(x_layer, y_layer, 1); 1255 1256 // check and get above maximum: 1257 float delta_x_above = 0, delta_y_above = 0; 1258 float max_above = getScoreMaxAbove(layer, x_layer, y_layer, center, ismax, delta_x_above, delta_y_above); 1259 1260 if (!ismax) 1261 return 0.0f; 1262 1263 float max; // to be returned 1264 1265 if (layer % 2 == 0) 1266 { // on octave 1267 // treat the patch below: 1268 float delta_x_below, delta_y_below; 1269 float max_below_float; 1270 int max_below = 0; 1271 if (layer == 0) 1272 { 1273 // guess the lower intra octave... 1274 const BriskLayer& l = pyramid_[0]; 1275 int s_0_0 = l.getAgastScore_5_8(x_layer - 1, y_layer - 1, 1); 1276 max_below = s_0_0; 1277 int s_1_0 = l.getAgastScore_5_8(x_layer, y_layer - 1, 1); 1278 max_below = std::max(s_1_0, max_below); 1279 int s_2_0 = l.getAgastScore_5_8(x_layer + 1, y_layer - 1, 1); 1280 max_below = std::max(s_2_0, max_below); 1281 int s_2_1 = l.getAgastScore_5_8(x_layer + 1, y_layer, 1); 1282 max_below = std::max(s_2_1, max_below); 1283 int s_1_1 = l.getAgastScore_5_8(x_layer, y_layer, 1); 1284 max_below = std::max(s_1_1, max_below); 1285 int s_0_1 = l.getAgastScore_5_8(x_layer - 1, y_layer, 1); 1286 max_below = std::max(s_0_1, max_below); 1287 int s_0_2 = l.getAgastScore_5_8(x_layer - 1, y_layer + 1, 1); 1288 max_below = std::max(s_0_2, max_below); 1289 int s_1_2 = l.getAgastScore_5_8(x_layer, y_layer + 1, 1); 1290 max_below = std::max(s_1_2, max_below); 1291 int s_2_2 = l.getAgastScore_5_8(x_layer + 1, y_layer + 1, 1); 1292 max_below = std::max(s_2_2, max_below); 1293 1294 max_below_float = subpixel2D(s_0_0, s_0_1, s_0_2, s_1_0, s_1_1, s_1_2, s_2_0, s_2_1, s_2_2, delta_x_below, 1295 delta_y_below); 1296 max_below_float = (float)max_below; 1297 } 1298 else 1299 { 1300 max_below_float = getScoreMaxBelow(layer, x_layer, y_layer, center, ismax, delta_x_below, delta_y_below); 1301 if (!ismax) 1302 return 0; 1303 } 1304 1305 // get the patch on this layer: 1306 int s_0_0 = thisLayer.getAgastScore(x_layer - 1, y_layer - 1, 1); 1307 int s_1_0 = thisLayer.getAgastScore(x_layer, y_layer - 1, 1); 1308 int s_2_0 = thisLayer.getAgastScore(x_layer + 1, y_layer - 1, 1); 1309 int s_2_1 = thisLayer.getAgastScore(x_layer + 1, y_layer, 1); 1310 int s_1_1 = thisLayer.getAgastScore(x_layer, y_layer, 1); 1311 int s_0_1 = thisLayer.getAgastScore(x_layer - 1, y_layer, 1); 1312 int s_0_2 = thisLayer.getAgastScore(x_layer - 1, y_layer + 1, 1); 1313 int s_1_2 = thisLayer.getAgastScore(x_layer, y_layer + 1, 1); 1314 int s_2_2 = thisLayer.getAgastScore(x_layer + 1, y_layer + 1, 1); 1315 float delta_x_layer, delta_y_layer; 1316 float max_layer = subpixel2D(s_0_0, s_0_1, s_0_2, s_1_0, s_1_1, s_1_2, s_2_0, s_2_1, s_2_2, delta_x_layer, 1317 delta_y_layer); 1318 1319 // calculate the relative scale (1D maximum): 1320 if (layer == 0) 1321 { 1322 scale = refine1D_2(max_below_float, std::max(float(center), max_layer), max_above, max); 1323 } 1324 else 1325 scale = refine1D(max_below_float, std::max(float(center), max_layer), max_above, max); 1326 1327 if (scale > 1.0) 1328 { 1329 // interpolate the position: 1330 const float r0 = (1.5f - scale) / .5f; 1331 const float r1 = 1.0f - r0; 1332 x = (r0 * delta_x_layer + r1 * delta_x_above + float(x_layer)) * thisLayer.scale() + thisLayer.offset(); 1333 y = (r0 * delta_y_layer + r1 * delta_y_above + float(y_layer)) * thisLayer.scale() + thisLayer.offset(); 1334 } 1335 else 1336 { 1337 if (layer == 0) 1338 { 1339 // interpolate the position: 1340 const float r0 = (scale - 0.5f) / 0.5f; 1341 const float r_1 = 1.0f - r0; 1342 x = r0 * delta_x_layer + r_1 * delta_x_below + float(x_layer); 1343 y = r0 * delta_y_layer + r_1 * delta_y_below + float(y_layer); 1344 } 1345 else 1346 { 1347 // interpolate the position: 1348 const float r0 = (scale - 0.75f) / 0.25f; 1349 const float r_1 = 1.0f - r0; 1350 x = (r0 * delta_x_layer + r_1 * delta_x_below + float(x_layer)) * thisLayer.scale() + thisLayer.offset(); 1351 y = (r0 * delta_y_layer + r_1 * delta_y_below + float(y_layer)) * thisLayer.scale() + thisLayer.offset(); 1352 } 1353 } 1354 } 1355 else 1356 { 1357 // on intra 1358 // check the patch below: 1359 float delta_x_below, delta_y_below; 1360 float max_below = getScoreMaxBelow(layer, x_layer, y_layer, center, ismax, delta_x_below, delta_y_below); 1361 if (!ismax) 1362 return 0.0f; 1363 1364 // get the patch on this layer: 1365 int s_0_0 = thisLayer.getAgastScore(x_layer - 1, y_layer - 1, 1); 1366 int s_1_0 = thisLayer.getAgastScore(x_layer, y_layer - 1, 1); 1367 int s_2_0 = thisLayer.getAgastScore(x_layer + 1, y_layer - 1, 1); 1368 int s_2_1 = thisLayer.getAgastScore(x_layer + 1, y_layer, 1); 1369 int s_1_1 = thisLayer.getAgastScore(x_layer, y_layer, 1); 1370 int s_0_1 = thisLayer.getAgastScore(x_layer - 1, y_layer, 1); 1371 int s_0_2 = thisLayer.getAgastScore(x_layer - 1, y_layer + 1, 1); 1372 int s_1_2 = thisLayer.getAgastScore(x_layer, y_layer + 1, 1); 1373 int s_2_2 = thisLayer.getAgastScore(x_layer + 1, y_layer + 1, 1); 1374 float delta_x_layer, delta_y_layer; 1375 float max_layer = subpixel2D(s_0_0, s_0_1, s_0_2, s_1_0, s_1_1, s_1_2, s_2_0, s_2_1, s_2_2, delta_x_layer, 1376 delta_y_layer); 1377 1378 // calculate the relative scale (1D maximum): 1379 scale = refine1D_1(max_below, std::max(float(center), max_layer), max_above, max); 1380 if (scale > 1.0) 1381 { 1382 // interpolate the position: 1383 const float r0 = 4.0f - scale * 3.0f; 1384 const float r1 = 1.0f - r0; 1385 x = (r0 * delta_x_layer + r1 * delta_x_above + float(x_layer)) * thisLayer.scale() + thisLayer.offset(); 1386 y = (r0 * delta_y_layer + r1 * delta_y_above + float(y_layer)) * thisLayer.scale() + thisLayer.offset(); 1387 } 1388 else 1389 { 1390 // interpolate the position: 1391 const float r0 = scale * 3.0f - 2.0f; 1392 const float r_1 = 1.0f - r0; 1393 x = (r0 * delta_x_layer + r_1 * delta_x_below + float(x_layer)) * thisLayer.scale() + thisLayer.offset(); 1394 y = (r0 * delta_y_layer + r_1 * delta_y_below + float(y_layer)) * thisLayer.scale() + thisLayer.offset(); 1395 } 1396 } 1397 1398 // calculate the absolute scale: 1399 scale *= thisLayer.scale(); 1400 1401 // that's it, return the refined maximum: 1402 return max; 1403 } 1404 1405 // return the maximum of score patches above or below 1406 inline float 1407 BriskScaleSpace::getScoreMaxAbove(const int layer, const int x_layer, const int y_layer, const int threshold, 1408 bool& ismax, float& dx, float& dy) const 1409 { 1410 1411 ismax = false; 1412 // relevant floating point coordinates 1413 float x_1; 1414 float x1; 1415 float y_1; 1416 float y1; 1417 1418 // the layer above 1419 CV_Assert(layer + 1 < layers_); 1420 const BriskLayer& layerAbove = pyramid_[layer + 1]; 1421 1422 if (layer % 2 == 0) 1423 { 1424 // octave 1425 x_1 = float(4 * (x_layer) - 1 - 2) / 6.0f; 1426 x1 = float(4 * (x_layer) - 1 + 2) / 6.0f; 1427 y_1 = float(4 * (y_layer) - 1 - 2) / 6.0f; 1428 y1 = float(4 * (y_layer) - 1 + 2) / 6.0f; 1429 } 1430 else 1431 { 1432 // intra 1433 x_1 = float(6 * (x_layer) - 1 - 3) / 8.0f; 1434 x1 = float(6 * (x_layer) - 1 + 3) / 8.0f; 1435 y_1 = float(6 * (y_layer) - 1 - 3) / 8.0f; 1436 y1 = float(6 * (y_layer) - 1 + 3) / 8.0f; 1437 } 1438 1439 // check the first row 1440 int max_x = (int)x_1 + 1; 1441 int max_y = (int)y_1 + 1; 1442 float tmp_max; 1443 float maxval = (float)layerAbove.getAgastScore(x_1, y_1, 1); 1444 if (maxval > threshold) 1445 return 0; 1446 for (int x = (int)x_1 + 1; x <= int(x1); x++) 1447 { 1448 tmp_max = (float)layerAbove.getAgastScore(float(x), y_1, 1); 1449 if (tmp_max > threshold) 1450 return 0; 1451 if (tmp_max > maxval) 1452 { 1453 maxval = tmp_max; 1454 max_x = x; 1455 } 1456 } 1457 tmp_max = (float)layerAbove.getAgastScore(x1, y_1, 1); 1458 if (tmp_max > threshold) 1459 return 0; 1460 if (tmp_max > maxval) 1461 { 1462 maxval = tmp_max; 1463 max_x = int(x1); 1464 } 1465 1466 // middle rows 1467 for (int y = (int)y_1 + 1; y <= int(y1); y++) 1468 { 1469 tmp_max = (float)layerAbove.getAgastScore(x_1, float(y), 1); 1470 if (tmp_max > threshold) 1471 return 0; 1472 if (tmp_max > maxval) 1473 { 1474 maxval = tmp_max; 1475 max_x = int(x_1 + 1); 1476 max_y = y; 1477 } 1478 for (int x = (int)x_1 + 1; x <= int(x1); x++) 1479 { 1480 tmp_max = (float)layerAbove.getAgastScore(x, y, 1); 1481 if (tmp_max > threshold) 1482 return 0; 1483 if (tmp_max > maxval) 1484 { 1485 maxval = tmp_max; 1486 max_x = x; 1487 max_y = y; 1488 } 1489 } 1490 tmp_max = (float)layerAbove.getAgastScore(x1, float(y), 1); 1491 if (tmp_max > threshold) 1492 return 0; 1493 if (tmp_max > maxval) 1494 { 1495 maxval = tmp_max; 1496 max_x = int(x1); 1497 max_y = y; 1498 } 1499 } 1500 1501 // bottom row 1502 tmp_max = (float)layerAbove.getAgastScore(x_1, y1, 1); 1503 if (tmp_max > maxval) 1504 { 1505 maxval = tmp_max; 1506 max_x = int(x_1 + 1); 1507 max_y = int(y1); 1508 } 1509 for (int x = (int)x_1 + 1; x <= int(x1); x++) 1510 { 1511 tmp_max = (float)layerAbove.getAgastScore(float(x), y1, 1); 1512 if (tmp_max > maxval) 1513 { 1514 maxval = tmp_max; 1515 max_x = x; 1516 max_y = int(y1); 1517 } 1518 } 1519 tmp_max = (float)layerAbove.getAgastScore(x1, y1, 1); 1520 if (tmp_max > maxval) 1521 { 1522 maxval = tmp_max; 1523 max_x = int(x1); 1524 max_y = int(y1); 1525 } 1526 1527 //find dx/dy: 1528 int s_0_0 = layerAbove.getAgastScore(max_x - 1, max_y - 1, 1); 1529 int s_1_0 = layerAbove.getAgastScore(max_x, max_y - 1, 1); 1530 int s_2_0 = layerAbove.getAgastScore(max_x + 1, max_y - 1, 1); 1531 int s_2_1 = layerAbove.getAgastScore(max_x + 1, max_y, 1); 1532 int s_1_1 = layerAbove.getAgastScore(max_x, max_y, 1); 1533 int s_0_1 = layerAbove.getAgastScore(max_x - 1, max_y, 1); 1534 int s_0_2 = layerAbove.getAgastScore(max_x - 1, max_y + 1, 1); 1535 int s_1_2 = layerAbove.getAgastScore(max_x, max_y + 1, 1); 1536 int s_2_2 = layerAbove.getAgastScore(max_x + 1, max_y + 1, 1); 1537 float dx_1, dy_1; 1538 float refined_max = subpixel2D(s_0_0, s_0_1, s_0_2, s_1_0, s_1_1, s_1_2, s_2_0, s_2_1, s_2_2, dx_1, dy_1); 1539 1540 // calculate dx/dy in above coordinates 1541 float real_x = float(max_x) + dx_1; 1542 float real_y = float(max_y) + dy_1; 1543 bool returnrefined = true; 1544 if (layer % 2 == 0) 1545 { 1546 dx = (real_x * 6.0f + 1.0f) / 4.0f - float(x_layer); 1547 dy = (real_y * 6.0f + 1.0f) / 4.0f - float(y_layer); 1548 } 1549 else 1550 { 1551 dx = (real_x * 8.0f + 1.0f) / 6.0f - float(x_layer); 1552 dy = (real_y * 8.0f + 1.0f) / 6.0f - float(y_layer); 1553 } 1554 1555 // saturate 1556 if (dx > 1.0f) 1557 { 1558 dx = 1.0f; 1559 returnrefined = false; 1560 } 1561 if (dx < -1.0f) 1562 { 1563 dx = -1.0f; 1564 returnrefined = false; 1565 } 1566 if (dy > 1.0f) 1567 { 1568 dy = 1.0f; 1569 returnrefined = false; 1570 } 1571 if (dy < -1.0f) 1572 { 1573 dy = -1.0f; 1574 returnrefined = false; 1575 } 1576 1577 // done and ok. 1578 ismax = true; 1579 if (returnrefined) 1580 { 1581 return std::max(refined_max, maxval); 1582 } 1583 return maxval; 1584 } 1585 1586 inline float 1587 BriskScaleSpace::getScoreMaxBelow(const int layer, const int x_layer, const int y_layer, const int threshold, 1588 bool& ismax, float& dx, float& dy) const 1589 { 1590 ismax = false; 1591 1592 // relevant floating point coordinates 1593 float x_1; 1594 float x1; 1595 float y_1; 1596 float y1; 1597 1598 if (layer % 2 == 0) 1599 { 1600 // octave 1601 x_1 = float(8 * (x_layer) + 1 - 4) / 6.0f; 1602 x1 = float(8 * (x_layer) + 1 + 4) / 6.0f; 1603 y_1 = float(8 * (y_layer) + 1 - 4) / 6.0f; 1604 y1 = float(8 * (y_layer) + 1 + 4) / 6.0f; 1605 } 1606 else 1607 { 1608 x_1 = float(6 * (x_layer) + 1 - 3) / 4.0f; 1609 x1 = float(6 * (x_layer) + 1 + 3) / 4.0f; 1610 y_1 = float(6 * (y_layer) + 1 - 3) / 4.0f; 1611 y1 = float(6 * (y_layer) + 1 + 3) / 4.0f; 1612 } 1613 1614 // the layer below 1615 CV_Assert(layer > 0); 1616 const BriskLayer& layerBelow = pyramid_[layer - 1]; 1617 1618 // check the first row 1619 int max_x = (int)x_1 + 1; 1620 int max_y = (int)y_1 + 1; 1621 float tmp_max; 1622 float max = (float)layerBelow.getAgastScore(x_1, y_1, 1); 1623 if (max > threshold) 1624 return 0; 1625 for (int x = (int)x_1 + 1; x <= int(x1); x++) 1626 { 1627 tmp_max = (float)layerBelow.getAgastScore(float(x), y_1, 1); 1628 if (tmp_max > threshold) 1629 return 0; 1630 if (tmp_max > max) 1631 { 1632 max = tmp_max; 1633 max_x = x; 1634 } 1635 } 1636 tmp_max = (float)layerBelow.getAgastScore(x1, y_1, 1); 1637 if (tmp_max > threshold) 1638 return 0; 1639 if (tmp_max > max) 1640 { 1641 max = tmp_max; 1642 max_x = int(x1); 1643 } 1644 1645 // middle rows 1646 for (int y = (int)y_1 + 1; y <= int(y1); y++) 1647 { 1648 tmp_max = (float)layerBelow.getAgastScore(x_1, float(y), 1); 1649 if (tmp_max > threshold) 1650 return 0; 1651 if (tmp_max > max) 1652 { 1653 max = tmp_max; 1654 max_x = int(x_1 + 1); 1655 max_y = y; 1656 } 1657 for (int x = (int)x_1 + 1; x <= int(x1); x++) 1658 { 1659 tmp_max = (float)layerBelow.getAgastScore(x, y, 1); 1660 if (tmp_max > threshold) 1661 return 0; 1662 if (tmp_max == max) 1663 { 1664 const int t1 = 2 1665 * (layerBelow.getAgastScore(x - 1, y, 1) + layerBelow.getAgastScore(x + 1, y, 1) 1666 + layerBelow.getAgastScore(x, y + 1, 1) + layerBelow.getAgastScore(x, y - 1, 1)) 1667 + (layerBelow.getAgastScore(x + 1, y + 1, 1) + layerBelow.getAgastScore(x - 1, y + 1, 1) 1668 + layerBelow.getAgastScore(x + 1, y - 1, 1) + layerBelow.getAgastScore(x - 1, y - 1, 1)); 1669 const int t2 = 2 1670 * (layerBelow.getAgastScore(max_x - 1, max_y, 1) + layerBelow.getAgastScore(max_x + 1, max_y, 1) 1671 + layerBelow.getAgastScore(max_x, max_y + 1, 1) + layerBelow.getAgastScore(max_x, max_y - 1, 1)) 1672 + (layerBelow.getAgastScore(max_x + 1, max_y + 1, 1) + layerBelow.getAgastScore(max_x - 1, 1673 max_y + 1, 1) 1674 + layerBelow.getAgastScore(max_x + 1, max_y - 1, 1) 1675 + layerBelow.getAgastScore(max_x - 1, max_y - 1, 1)); 1676 if (t1 > t2) 1677 { 1678 max_x = x; 1679 max_y = y; 1680 } 1681 } 1682 if (tmp_max > max) 1683 { 1684 max = tmp_max; 1685 max_x = x; 1686 max_y = y; 1687 } 1688 } 1689 tmp_max = (float)layerBelow.getAgastScore(x1, float(y), 1); 1690 if (tmp_max > threshold) 1691 return 0; 1692 if (tmp_max > max) 1693 { 1694 max = tmp_max; 1695 max_x = int(x1); 1696 max_y = y; 1697 } 1698 } 1699 1700 // bottom row 1701 tmp_max = (float)layerBelow.getAgastScore(x_1, y1, 1); 1702 if (tmp_max > max) 1703 { 1704 max = tmp_max; 1705 max_x = int(x_1 + 1); 1706 max_y = int(y1); 1707 } 1708 for (int x = (int)x_1 + 1; x <= int(x1); x++) 1709 { 1710 tmp_max = (float)layerBelow.getAgastScore(float(x), y1, 1); 1711 if (tmp_max > max) 1712 { 1713 max = tmp_max; 1714 max_x = x; 1715 max_y = int(y1); 1716 } 1717 } 1718 tmp_max = (float)layerBelow.getAgastScore(x1, y1, 1); 1719 if (tmp_max > max) 1720 { 1721 max = tmp_max; 1722 max_x = int(x1); 1723 max_y = int(y1); 1724 } 1725 1726 //find dx/dy: 1727 int s_0_0 = layerBelow.getAgastScore(max_x - 1, max_y - 1, 1); 1728 int s_1_0 = layerBelow.getAgastScore(max_x, max_y - 1, 1); 1729 int s_2_0 = layerBelow.getAgastScore(max_x + 1, max_y - 1, 1); 1730 int s_2_1 = layerBelow.getAgastScore(max_x + 1, max_y, 1); 1731 int s_1_1 = layerBelow.getAgastScore(max_x, max_y, 1); 1732 int s_0_1 = layerBelow.getAgastScore(max_x - 1, max_y, 1); 1733 int s_0_2 = layerBelow.getAgastScore(max_x - 1, max_y + 1, 1); 1734 int s_1_2 = layerBelow.getAgastScore(max_x, max_y + 1, 1); 1735 int s_2_2 = layerBelow.getAgastScore(max_x + 1, max_y + 1, 1); 1736 float dx_1, dy_1; 1737 float refined_max = subpixel2D(s_0_0, s_0_1, s_0_2, s_1_0, s_1_1, s_1_2, s_2_0, s_2_1, s_2_2, dx_1, dy_1); 1738 1739 // calculate dx/dy in above coordinates 1740 float real_x = float(max_x) + dx_1; 1741 float real_y = float(max_y) + dy_1; 1742 bool returnrefined = true; 1743 if (layer % 2 == 0) 1744 { 1745 dx = (float)((real_x * 6.0 + 1.0) / 8.0) - float(x_layer); 1746 dy = (float)((real_y * 6.0 + 1.0) / 8.0) - float(y_layer); 1747 } 1748 else 1749 { 1750 dx = (float)((real_x * 4.0 - 1.0) / 6.0) - float(x_layer); 1751 dy = (float)((real_y * 4.0 - 1.0) / 6.0) - float(y_layer); 1752 } 1753 1754 // saturate 1755 if (dx > 1.0) 1756 { 1757 dx = 1.0f; 1758 returnrefined = false; 1759 } 1760 if (dx < -1.0f) 1761 { 1762 dx = -1.0f; 1763 returnrefined = false; 1764 } 1765 if (dy > 1.0f) 1766 { 1767 dy = 1.0f; 1768 returnrefined = false; 1769 } 1770 if (dy < -1.0f) 1771 { 1772 dy = -1.0f; 1773 returnrefined = false; 1774 } 1775 1776 // done and ok. 1777 ismax = true; 1778 if (returnrefined) 1779 { 1780 return std::max(refined_max, max); 1781 } 1782 return max; 1783 } 1784 1785 inline float 1786 BriskScaleSpace::refine1D(const float s_05, const float s0, const float s05, float& max) const 1787 { 1788 int i_05 = int(1024.0 * s_05 + 0.5); 1789 int i0 = int(1024.0 * s0 + 0.5); 1790 int i05 = int(1024.0 * s05 + 0.5); 1791 1792 // 16.0000 -24.0000 8.0000 1793 // -40.0000 54.0000 -14.0000 1794 // 24.0000 -27.0000 6.0000 1795 1796 int three_a = 16 * i_05 - 24 * i0 + 8 * i05; 1797 // second derivative must be negative: 1798 if (three_a >= 0) 1799 { 1800 if (s0 >= s_05 && s0 >= s05) 1801 { 1802 max = s0; 1803 return 1.0f; 1804 } 1805 if (s_05 >= s0 && s_05 >= s05) 1806 { 1807 max = s_05; 1808 return 0.75f; 1809 } 1810 if (s05 >= s0 && s05 >= s_05) 1811 { 1812 max = s05; 1813 return 1.5f; 1814 } 1815 } 1816 1817 int three_b = -40 * i_05 + 54 * i0 - 14 * i05; 1818 // calculate max location: 1819 float ret_val = -float(three_b) / float(2 * three_a); 1820 // saturate and return 1821 if (ret_val < 0.75) 1822 ret_val = 0.75; 1823 else if (ret_val > 1.5) 1824 ret_val = 1.5; // allow to be slightly off bounds ...? 1825 int three_c = +24 * i_05 - 27 * i0 + 6 * i05; 1826 max = float(three_c) + float(three_a) * ret_val * ret_val + float(three_b) * ret_val; 1827 max /= 3072.0f; 1828 return ret_val; 1829 } 1830 1831 inline float 1832 BriskScaleSpace::refine1D_1(const float s_05, const float s0, const float s05, float& max) const 1833 { 1834 int i_05 = int(1024.0 * s_05 + 0.5); 1835 int i0 = int(1024.0 * s0 + 0.5); 1836 int i05 = int(1024.0 * s05 + 0.5); 1837 1838 // 4.5000 -9.0000 4.5000 1839 //-10.5000 18.0000 -7.5000 1840 // 6.0000 -8.0000 3.0000 1841 1842 int two_a = 9 * i_05 - 18 * i0 + 9 * i05; 1843 // second derivative must be negative: 1844 if (two_a >= 0) 1845 { 1846 if (s0 >= s_05 && s0 >= s05) 1847 { 1848 max = s0; 1849 return 1.0f; 1850 } 1851 if (s_05 >= s0 && s_05 >= s05) 1852 { 1853 max = s_05; 1854 return 0.6666666666666666666666666667f; 1855 } 1856 if (s05 >= s0 && s05 >= s_05) 1857 { 1858 max = s05; 1859 return 1.3333333333333333333333333333f; 1860 } 1861 } 1862 1863 int two_b = -21 * i_05 + 36 * i0 - 15 * i05; 1864 // calculate max location: 1865 float ret_val = -float(two_b) / float(2 * two_a); 1866 // saturate and return 1867 if (ret_val < 0.6666666666666666666666666667f) 1868 ret_val = 0.666666666666666666666666667f; 1869 else if (ret_val > 1.33333333333333333333333333f) 1870 ret_val = 1.333333333333333333333333333f; 1871 int two_c = +12 * i_05 - 16 * i0 + 6 * i05; 1872 max = float(two_c) + float(two_a) * ret_val * ret_val + float(two_b) * ret_val; 1873 max /= 2048.0f; 1874 return ret_val; 1875 } 1876 1877 inline float 1878 BriskScaleSpace::refine1D_2(const float s_05, const float s0, const float s05, float& max) const 1879 { 1880 int i_05 = int(1024.0 * s_05 + 0.5); 1881 int i0 = int(1024.0 * s0 + 0.5); 1882 int i05 = int(1024.0 * s05 + 0.5); 1883 1884 // 18.0000 -30.0000 12.0000 1885 // -45.0000 65.0000 -20.0000 1886 // 27.0000 -30.0000 8.0000 1887 1888 int a = 2 * i_05 - 4 * i0 + 2 * i05; 1889 // second derivative must be negative: 1890 if (a >= 0) 1891 { 1892 if (s0 >= s_05 && s0 >= s05) 1893 { 1894 max = s0; 1895 return 1.0f; 1896 } 1897 if (s_05 >= s0 && s_05 >= s05) 1898 { 1899 max = s_05; 1900 return 0.7f; 1901 } 1902 if (s05 >= s0 && s05 >= s_05) 1903 { 1904 max = s05; 1905 return 1.5f; 1906 } 1907 } 1908 1909 int b = -5 * i_05 + 8 * i0 - 3 * i05; 1910 // calculate max location: 1911 float ret_val = -float(b) / float(2 * a); 1912 // saturate and return 1913 if (ret_val < 0.7f) 1914 ret_val = 0.7f; 1915 else if (ret_val > 1.5f) 1916 ret_val = 1.5f; // allow to be slightly off bounds ...? 1917 int c = +3 * i_05 - 3 * i0 + 1 * i05; 1918 max = float(c) + float(a) * ret_val * ret_val + float(b) * ret_val; 1919 max /= 1024; 1920 return ret_val; 1921 } 1922 1923 inline float 1924 BriskScaleSpace::subpixel2D(const int s_0_0, const int s_0_1, const int s_0_2, const int s_1_0, const int s_1_1, 1925 const int s_1_2, const int s_2_0, const int s_2_1, const int s_2_2, float& delta_x, 1926 float& delta_y) const 1927 { 1928 1929 // the coefficients of the 2d quadratic function least-squares fit: 1930 int tmp1 = s_0_0 + s_0_2 - 2 * s_1_1 + s_2_0 + s_2_2; 1931 int coeff1 = 3 * (tmp1 + s_0_1 - ((s_1_0 + s_1_2) << 1) + s_2_1); 1932 int coeff2 = 3 * (tmp1 - ((s_0_1 + s_2_1) << 1) + s_1_0 + s_1_2); 1933 int tmp2 = s_0_2 - s_2_0; 1934 int tmp3 = (s_0_0 + tmp2 - s_2_2); 1935 int tmp4 = tmp3 - 2 * tmp2; 1936 int coeff3 = -3 * (tmp3 + s_0_1 - s_2_1); 1937 int coeff4 = -3 * (tmp4 + s_1_0 - s_1_2); 1938 int coeff5 = (s_0_0 - s_0_2 - s_2_0 + s_2_2) << 2; 1939 int coeff6 = -(s_0_0 + s_0_2 - ((s_1_0 + s_0_1 + s_1_2 + s_2_1) << 1) - 5 * s_1_1 + s_2_0 + s_2_2) << 1; 1940 1941 // 2nd derivative test: 1942 int H_det = 4 * coeff1 * coeff2 - coeff5 * coeff5; 1943 1944 if (H_det == 0) 1945 { 1946 delta_x = 0.0f; 1947 delta_y = 0.0f; 1948 return float(coeff6) / 18.0f; 1949 } 1950 1951 if (!(H_det > 0 && coeff1 < 0)) 1952 { 1953 // The maximum must be at the one of the 4 patch corners. 1954 int tmp_max = coeff3 + coeff4 + coeff5; 1955 delta_x = 1.0f; 1956 delta_y = 1.0f; 1957 1958 int tmp = -coeff3 + coeff4 - coeff5; 1959 if (tmp > tmp_max) 1960 { 1961 tmp_max = tmp; 1962 delta_x = -1.0f; 1963 delta_y = 1.0f; 1964 } 1965 tmp = coeff3 - coeff4 - coeff5; 1966 if (tmp > tmp_max) 1967 { 1968 tmp_max = tmp; 1969 delta_x = 1.0f; 1970 delta_y = -1.0f; 1971 } 1972 tmp = -coeff3 - coeff4 + coeff5; 1973 if (tmp > tmp_max) 1974 { 1975 tmp_max = tmp; 1976 delta_x = -1.0f; 1977 delta_y = -1.0f; 1978 } 1979 return float(tmp_max + coeff1 + coeff2 + coeff6) / 18.0f; 1980 } 1981 1982 // this is hopefully the normal outcome of the Hessian test 1983 delta_x = float(2 * coeff2 * coeff3 - coeff4 * coeff5) / float(-H_det); 1984 delta_y = float(2 * coeff1 * coeff4 - coeff3 * coeff5) / float(-H_det); 1985 // TODO: this is not correct, but easy, so perform a real boundary maximum search: 1986 bool tx = false; 1987 bool tx_ = false; 1988 bool ty = false; 1989 bool ty_ = false; 1990 if (delta_x > 1.0) 1991 tx = true; 1992 else if (delta_x < -1.0) 1993 tx_ = true; 1994 if (delta_y > 1.0) 1995 ty = true; 1996 if (delta_y < -1.0) 1997 ty_ = true; 1998 1999 if (tx || tx_ || ty || ty_) 2000 { 2001 // get two candidates: 2002 float delta_x1 = 0.0f, delta_x2 = 0.0f, delta_y1 = 0.0f, delta_y2 = 0.0f; 2003 if (tx) 2004 { 2005 delta_x1 = 1.0f; 2006 delta_y1 = -float(coeff4 + coeff5) / float(2 * coeff2); 2007 if (delta_y1 > 1.0f) 2008 delta_y1 = 1.0f; 2009 else if (delta_y1 < -1.0f) 2010 delta_y1 = -1.0f; 2011 } 2012 else if (tx_) 2013 { 2014 delta_x1 = -1.0f; 2015 delta_y1 = -float(coeff4 - coeff5) / float(2 * coeff2); 2016 if (delta_y1 > 1.0f) 2017 delta_y1 = 1.0f; 2018 else if (delta_y1 < -1.0) 2019 delta_y1 = -1.0f; 2020 } 2021 if (ty) 2022 { 2023 delta_y2 = 1.0f; 2024 delta_x2 = -float(coeff3 + coeff5) / float(2 * coeff1); 2025 if (delta_x2 > 1.0f) 2026 delta_x2 = 1.0f; 2027 else if (delta_x2 < -1.0f) 2028 delta_x2 = -1.0f; 2029 } 2030 else if (ty_) 2031 { 2032 delta_y2 = -1.0f; 2033 delta_x2 = -float(coeff3 - coeff5) / float(2 * coeff1); 2034 if (delta_x2 > 1.0f) 2035 delta_x2 = 1.0f; 2036 else if (delta_x2 < -1.0f) 2037 delta_x2 = -1.0f; 2038 } 2039 // insert both options for evaluation which to pick 2040 float max1 = (coeff1 * delta_x1 * delta_x1 + coeff2 * delta_y1 * delta_y1 + coeff3 * delta_x1 + coeff4 * delta_y1 2041 + coeff5 * delta_x1 * delta_y1 + coeff6) 2042 / 18.0f; 2043 float max2 = (coeff1 * delta_x2 * delta_x2 + coeff2 * delta_y2 * delta_y2 + coeff3 * delta_x2 + coeff4 * delta_y2 2044 + coeff5 * delta_x2 * delta_y2 + coeff6) 2045 / 18.0f; 2046 if (max1 > max2) 2047 { 2048 delta_x = delta_x1; 2049 delta_y = delta_x1; 2050 return max1; 2051 } 2052 else 2053 { 2054 delta_x = delta_x2; 2055 delta_y = delta_x2; 2056 return max2; 2057 } 2058 } 2059 2060 // this is the case of the maximum inside the boundaries: 2061 return (coeff1 * delta_x * delta_x + coeff2 * delta_y * delta_y + coeff3 * delta_x + coeff4 * delta_y 2062 + coeff5 * delta_x * delta_y + coeff6) 2063 / 18.0f; 2064 } 2065 2066 // construct a layer 2067 BriskLayer::BriskLayer(const cv::Mat& img_in, float scale_in, float offset_in) 2068 { 2069 img_ = img_in; 2070 scores_ = cv::Mat_<uchar>::zeros(img_in.rows, img_in.cols); 2071 // attention: this means that the passed image reference must point to persistent memory 2072 scale_ = scale_in; 2073 offset_ = offset_in; 2074 // create an agast detector 2075 oast_9_16_ = AgastFeatureDetector::create(1, false, AgastFeatureDetector::OAST_9_16); 2076 makeAgastOffsets(pixel_5_8_, (int)img_.step, AgastFeatureDetector::AGAST_5_8); 2077 makeAgastOffsets(pixel_9_16_, (int)img_.step, AgastFeatureDetector::OAST_9_16); 2078 } 2079 // derive a layer 2080 BriskLayer::BriskLayer(const BriskLayer& layer, int mode) 2081 { 2082 if (mode == CommonParams::HALFSAMPLE) 2083 { 2084 img_.create(layer.img().rows / 2, layer.img().cols / 2, CV_8U); 2085 halfsample(layer.img(), img_); 2086 scale_ = layer.scale() * 2; 2087 offset_ = 0.5f * scale_ - 0.5f; 2088 } 2089 else 2090 { 2091 img_.create(2 * (layer.img().rows / 3), 2 * (layer.img().cols / 3), CV_8U); 2092 twothirdsample(layer.img(), img_); 2093 scale_ = layer.scale() * 1.5f; 2094 offset_ = 0.5f * scale_ - 0.5f; 2095 } 2096 scores_ = cv::Mat::zeros(img_.rows, img_.cols, CV_8U); 2097 oast_9_16_ = AgastFeatureDetector::create(1, false, AgastFeatureDetector::OAST_9_16); 2098 makeAgastOffsets(pixel_5_8_, (int)img_.step, AgastFeatureDetector::AGAST_5_8); 2099 makeAgastOffsets(pixel_9_16_, (int)img_.step, AgastFeatureDetector::OAST_9_16); 2100 } 2101 2102 // Agast 2103 // wraps the agast class 2104 void 2105 BriskLayer::getAgastPoints(int threshold, std::vector<KeyPoint>& keypoints) 2106 { 2107 oast_9_16_->setThreshold(threshold); 2108 oast_9_16_->detect(img_, keypoints); 2109 2110 // also write scores 2111 const size_t num = keypoints.size(); 2112 2113 for (size_t i = 0; i < num; i++) 2114 scores_((int)keypoints[i].pt.y, (int)keypoints[i].pt.x) = saturate_cast<uchar>(keypoints[i].response); 2115 } 2116 2117 inline int 2118 BriskLayer::getAgastScore(int x, int y, int threshold) const 2119 { 2120 if (x < 3 || y < 3) 2121 return 0; 2122 if (x >= img_.cols - 3 || y >= img_.rows - 3) 2123 return 0; 2124 uchar& score = (uchar&)scores_(y, x); 2125 if (score > 2) 2126 { 2127 return score; 2128 } 2129 score = (uchar)agast_cornerScore<AgastFeatureDetector::OAST_9_16>(&img_.at<uchar>(y, x), pixel_9_16_, threshold - 1); 2130 if (score < threshold) 2131 score = 0; 2132 return score; 2133 } 2134 2135 inline int 2136 BriskLayer::getAgastScore_5_8(int x, int y, int threshold) const 2137 { 2138 if (x < 2 || y < 2) 2139 return 0; 2140 if (x >= img_.cols - 2 || y >= img_.rows - 2) 2141 return 0; 2142 int score = agast_cornerScore<AgastFeatureDetector::AGAST_5_8>(&img_.at<uchar>(y, x), pixel_5_8_, threshold - 1); 2143 if (score < threshold) 2144 score = 0; 2145 return score; 2146 } 2147 2148 inline int 2149 BriskLayer::getAgastScore(float xf, float yf, int threshold_in, float scale_in) const 2150 { 2151 if (scale_in <= 1.0f) 2152 { 2153 // just do an interpolation inside the layer 2154 const int x = int(xf); 2155 const float rx1 = xf - float(x); 2156 const float rx = 1.0f - rx1; 2157 const int y = int(yf); 2158 const float ry1 = yf - float(y); 2159 const float ry = 1.0f - ry1; 2160 2161 return (uchar)(rx * ry * getAgastScore(x, y, threshold_in) + rx1 * ry * getAgastScore(x + 1, y, threshold_in) 2162 + rx * ry1 * getAgastScore(x, y + 1, threshold_in) + rx1 * ry1 * getAgastScore(x + 1, y + 1, threshold_in)); 2163 } 2164 else 2165 { 2166 // this means we overlap area smoothing 2167 const float halfscale = scale_in / 2.0f; 2168 // get the scores first: 2169 for (int x = int(xf - halfscale); x <= int(xf + halfscale + 1.0f); x++) 2170 { 2171 for (int y = int(yf - halfscale); y <= int(yf + halfscale + 1.0f); y++) 2172 { 2173 getAgastScore(x, y, threshold_in); 2174 } 2175 } 2176 // get the smoothed value 2177 return value(scores_, xf, yf, scale_in); 2178 } 2179 } 2180 2181 // access gray values (smoothed/interpolated) 2182 inline int 2183 BriskLayer::value(const cv::Mat& mat, float xf, float yf, float scale_in) const 2184 { 2185 CV_Assert(!mat.empty()); 2186 // get the position 2187 const int x = cvFloor(xf); 2188 const int y = cvFloor(yf); 2189 const cv::Mat& image = mat; 2190 const int& imagecols = image.cols; 2191 2192 // get the sigma_half: 2193 const float sigma_half = scale_in / 2; 2194 const float area = 4.0f * sigma_half * sigma_half; 2195 // calculate output: 2196 int ret_val; 2197 if (sigma_half < 0.5) 2198 { 2199 //interpolation multipliers: 2200 const int r_x = (int)((xf - x) * 1024); 2201 const int r_y = (int)((yf - y) * 1024); 2202 const int r_x_1 = (1024 - r_x); 2203 const int r_y_1 = (1024 - r_y); 2204 const uchar* ptr = image.ptr() + x + y * imagecols; 2205 // just interpolate: 2206 ret_val = (r_x_1 * r_y_1 * int(*ptr)); 2207 ptr++; 2208 ret_val += (r_x * r_y_1 * int(*ptr)); 2209 ptr += imagecols; 2210 ret_val += (r_x * r_y * int(*ptr)); 2211 ptr--; 2212 ret_val += (r_x_1 * r_y * int(*ptr)); 2213 return 0xFF & ((ret_val + 512) / 1024 / 1024); 2214 } 2215 2216 // this is the standard case (simple, not speed optimized yet): 2217 2218 // scaling: 2219 const int scaling = (int)(4194304.0f / area); 2220 const int scaling2 = (int)(float(scaling) * area / 1024.0f); 2221 2222 // calculate borders 2223 const float x_1 = xf - sigma_half; 2224 const float x1 = xf + sigma_half; 2225 const float y_1 = yf - sigma_half; 2226 const float y1 = yf + sigma_half; 2227 2228 const int x_left = int(x_1 + 0.5); 2229 const int y_top = int(y_1 + 0.5); 2230 const int x_right = int(x1 + 0.5); 2231 const int y_bottom = int(y1 + 0.5); 2232 2233 // overlap area - multiplication factors: 2234 const float r_x_1 = float(x_left) - x_1 + 0.5f; 2235 const float r_y_1 = float(y_top) - y_1 + 0.5f; 2236 const float r_x1 = x1 - float(x_right) + 0.5f; 2237 const float r_y1 = y1 - float(y_bottom) + 0.5f; 2238 const int dx = x_right - x_left - 1; 2239 const int dy = y_bottom - y_top - 1; 2240 const int A = (int)((r_x_1 * r_y_1) * scaling); 2241 const int B = (int)((r_x1 * r_y_1) * scaling); 2242 const int C = (int)((r_x1 * r_y1) * scaling); 2243 const int D = (int)((r_x_1 * r_y1) * scaling); 2244 const int r_x_1_i = (int)(r_x_1 * scaling); 2245 const int r_y_1_i = (int)(r_y_1 * scaling); 2246 const int r_x1_i = (int)(r_x1 * scaling); 2247 const int r_y1_i = (int)(r_y1 * scaling); 2248 2249 // now the calculation: 2250 const uchar* ptr = image.ptr() + x_left + imagecols * y_top; 2251 // first row: 2252 ret_val = A * int(*ptr); 2253 ptr++; 2254 const uchar* end1 = ptr + dx; 2255 for (; ptr < end1; ptr++) 2256 { 2257 ret_val += r_y_1_i * int(*ptr); 2258 } 2259 ret_val += B * int(*ptr); 2260 // middle ones: 2261 ptr += imagecols - dx - 1; 2262 const uchar* end_j = ptr + dy * imagecols; 2263 for (; ptr < end_j; ptr += imagecols - dx - 1) 2264 { 2265 ret_val += r_x_1_i * int(*ptr); 2266 ptr++; 2267 const uchar* end2 = ptr + dx; 2268 for (; ptr < end2; ptr++) 2269 { 2270 ret_val += int(*ptr) * scaling; 2271 } 2272 ret_val += r_x1_i * int(*ptr); 2273 } 2274 // last row: 2275 ret_val += D * int(*ptr); 2276 ptr++; 2277 const uchar* end3 = ptr + dx; 2278 for (; ptr < end3; ptr++) 2279 { 2280 ret_val += r_y1_i * int(*ptr); 2281 } 2282 ret_val += C * int(*ptr); 2283 2284 return 0xFF & ((ret_val + scaling2 / 2) / scaling2 / 1024); 2285 } 2286 2287 // half sampling 2288 inline void 2289 BriskLayer::halfsample(const cv::Mat& srcimg, cv::Mat& dstimg) 2290 { 2291 // make sure the destination image is of the right size: 2292 CV_Assert(srcimg.cols / 2 == dstimg.cols); 2293 CV_Assert(srcimg.rows / 2 == dstimg.rows); 2294 2295 // handle non-SSE case 2296 resize(srcimg, dstimg, dstimg.size(), 0, 0, INTER_AREA); 2297 } 2298 2299 inline void 2300 BriskLayer::twothirdsample(const cv::Mat& srcimg, cv::Mat& dstimg) 2301 { 2302 // make sure the destination image is of the right size: 2303 CV_Assert((srcimg.cols / 3) * 2 == dstimg.cols); 2304 CV_Assert((srcimg.rows / 3) * 2 == dstimg.rows); 2305 2306 resize(srcimg, dstimg, dstimg.size(), 0, 0, INTER_AREA); 2307 } 2308 2309 Ptr<BRISK> BRISK::create(int thresh, int octaves, float patternScale) 2310 { 2311 return makePtr<BRISK_Impl>(thresh, octaves, patternScale); 2312 } 2313 2314 // custom setup 2315 Ptr<BRISK> BRISK::create(const std::vector<float> &radiusList, const std::vector<int> &numberList, 2316 float dMax, float dMin, const std::vector<int>& indexChange) 2317 { 2318 return makePtr<BRISK_Impl>(radiusList, numberList, dMax, dMin, indexChange); 2319 } 2320 2321 } 2322
