1 /** 2 * @file AKAZEFeatures.cpp 3 * @brief Main class for detecting and describing binary features in an 4 * accelerated nonlinear scale space 5 * @date Sep 15, 2013 6 * @author Pablo F. Alcantarilla, Jesus Nuevo 7 */ 8 9 #include "../precomp.hpp" 10 #include "AKAZEFeatures.h" 11 #include "fed.h" 12 #include "nldiffusion_functions.h" 13 #include "utils.h" 14 15 #include <iostream> 16 17 // Namespaces 18 namespace cv 19 { 20 using namespace std; 21 22 /* ************************************************************************* */ 23 /** 24 * @brief AKAZEFeatures constructor with input options 25 * @param options AKAZEFeatures configuration options 26 * @note This constructor allocates memory for the nonlinear scale space 27 */ 28 AKAZEFeatures::AKAZEFeatures(const AKAZEOptions& options) : options_(options) { 29 30 ncycles_ = 0; 31 reordering_ = true; 32 33 if (options_.descriptor_size > 0 && options_.descriptor >= AKAZE::DESCRIPTOR_MLDB_UPRIGHT) { 34 generateDescriptorSubsample(descriptorSamples_, descriptorBits_, options_.descriptor_size, 35 options_.descriptor_pattern_size, options_.descriptor_channels); 36 } 37 38 Allocate_Memory_Evolution(); 39 } 40 41 /* ************************************************************************* */ 42 /** 43 * @brief This method allocates the memory for the nonlinear diffusion evolution 44 */ 45 void AKAZEFeatures::Allocate_Memory_Evolution(void) { 46 47 float rfactor = 0.0f; 48 int level_height = 0, level_width = 0; 49 50 // Allocate the dimension of the matrices for the evolution 51 for (int i = 0, power = 1; i <= options_.omax - 1; i++, power *= 2) { 52 rfactor = 1.0f / power; 53 level_height = (int)(options_.img_height*rfactor); 54 level_width = (int)(options_.img_width*rfactor); 55 56 // Smallest possible octave and allow one scale if the image is small 57 if ((level_width < 80 || level_height < 40) && i != 0) { 58 options_.omax = i; 59 break; 60 } 61 62 for (int j = 0; j < options_.nsublevels; j++) { 63 TEvolution step; 64 step.Lx = Mat::zeros(level_height, level_width, CV_32F); 65 step.Ly = Mat::zeros(level_height, level_width, CV_32F); 66 step.Lxx = Mat::zeros(level_height, level_width, CV_32F); 67 step.Lxy = Mat::zeros(level_height, level_width, CV_32F); 68 step.Lyy = Mat::zeros(level_height, level_width, CV_32F); 69 step.Lt = Mat::zeros(level_height, level_width, CV_32F); 70 step.Ldet = Mat::zeros(level_height, level_width, CV_32F); 71 step.Lsmooth = Mat::zeros(level_height, level_width, CV_32F); 72 step.esigma = options_.soffset*pow(2.f, (float)(j) / (float)(options_.nsublevels) + i); 73 step.sigma_size = fRound(step.esigma); 74 step.etime = 0.5f*(step.esigma*step.esigma); 75 step.octave = i; 76 step.sublevel = j; 77 evolution_.push_back(step); 78 } 79 } 80 81 // Allocate memory for the number of cycles and time steps 82 for (size_t i = 1; i < evolution_.size(); i++) { 83 int naux = 0; 84 vector<float> tau; 85 float ttime = 0.0f; 86 ttime = evolution_[i].etime - evolution_[i - 1].etime; 87 naux = fed_tau_by_process_time(ttime, 1, 0.25f, reordering_, tau); 88 nsteps_.push_back(naux); 89 tsteps_.push_back(tau); 90 ncycles_++; 91 } 92 } 93 94 /* ************************************************************************* */ 95 /** 96 * @brief This method creates the nonlinear scale space for a given image 97 * @param img Input image for which the nonlinear scale space needs to be created 98 * @return 0 if the nonlinear scale space was created successfully, -1 otherwise 99 */ 100 int AKAZEFeatures::Create_Nonlinear_Scale_Space(const Mat& img) 101 { 102 CV_Assert(evolution_.size() > 0); 103 104 // Copy the original image to the first level of the evolution 105 img.copyTo(evolution_[0].Lt); 106 gaussian_2D_convolution(evolution_[0].Lt, evolution_[0].Lt, 0, 0, options_.soffset); 107 evolution_[0].Lt.copyTo(evolution_[0].Lsmooth); 108 109 // Allocate memory for the flow and step images 110 Mat Lflow = Mat::zeros(evolution_[0].Lt.rows, evolution_[0].Lt.cols, CV_32F); 111 Mat Lstep = Mat::zeros(evolution_[0].Lt.rows, evolution_[0].Lt.cols, CV_32F); 112 113 // First compute the kcontrast factor 114 options_.kcontrast = compute_k_percentile(img, options_.kcontrast_percentile, 1.0f, options_.kcontrast_nbins, 0, 0); 115 116 // Now generate the rest of evolution levels 117 for (size_t i = 1; i < evolution_.size(); i++) { 118 119 if (evolution_[i].octave > evolution_[i - 1].octave) { 120 halfsample_image(evolution_[i - 1].Lt, evolution_[i].Lt); 121 options_.kcontrast = options_.kcontrast*0.75f; 122 123 // Allocate memory for the resized flow and step images 124 Lflow = Mat::zeros(evolution_[i].Lt.rows, evolution_[i].Lt.cols, CV_32F); 125 Lstep = Mat::zeros(evolution_[i].Lt.rows, evolution_[i].Lt.cols, CV_32F); 126 } 127 else { 128 evolution_[i - 1].Lt.copyTo(evolution_[i].Lt); 129 } 130 131 gaussian_2D_convolution(evolution_[i].Lt, evolution_[i].Lsmooth, 0, 0, 1.0f); 132 133 // Compute the Gaussian derivatives Lx and Ly 134 image_derivatives_scharr(evolution_[i].Lsmooth, evolution_[i].Lx, 1, 0); 135 image_derivatives_scharr(evolution_[i].Lsmooth, evolution_[i].Ly, 0, 1); 136 137 // Compute the conductivity equation 138 switch (options_.diffusivity) { 139 case KAZE::DIFF_PM_G1: 140 pm_g1(evolution_[i].Lx, evolution_[i].Ly, Lflow, options_.kcontrast); 141 break; 142 case KAZE::DIFF_PM_G2: 143 pm_g2(evolution_[i].Lx, evolution_[i].Ly, Lflow, options_.kcontrast); 144 break; 145 case KAZE::DIFF_WEICKERT: 146 weickert_diffusivity(evolution_[i].Lx, evolution_[i].Ly, Lflow, options_.kcontrast); 147 break; 148 case KAZE::DIFF_CHARBONNIER: 149 charbonnier_diffusivity(evolution_[i].Lx, evolution_[i].Ly, Lflow, options_.kcontrast); 150 break; 151 default: 152 CV_Error(options_.diffusivity, "Diffusivity is not supported"); 153 break; 154 } 155 156 // Perform FED n inner steps 157 for (int j = 0; j < nsteps_[i - 1]; j++) { 158 nld_step_scalar(evolution_[i].Lt, Lflow, Lstep, tsteps_[i - 1][j]); 159 } 160 } 161 162 return 0; 163 } 164 165 /* ************************************************************************* */ 166 /** 167 * @brief This method selects interesting keypoints through the nonlinear scale space 168 * @param kpts Vector of detected keypoints 169 */ 170 void AKAZEFeatures::Feature_Detection(std::vector<KeyPoint>& kpts) 171 { 172 kpts.clear(); 173 Compute_Determinant_Hessian_Response(); 174 Find_Scale_Space_Extrema(kpts); 175 Do_Subpixel_Refinement(kpts); 176 } 177 178 /* ************************************************************************* */ 179 class MultiscaleDerivativesAKAZEInvoker : public ParallelLoopBody 180 { 181 public: 182 explicit MultiscaleDerivativesAKAZEInvoker(std::vector<TEvolution>& ev, const AKAZEOptions& opt) 183 : evolution_(&ev) 184 , options_(opt) 185 { 186 } 187 188 void operator()(const Range& range) const 189 { 190 std::vector<TEvolution>& evolution = *evolution_; 191 192 for (int i = range.start; i < range.end; i++) 193 { 194 float ratio = (float)fastpow(2, evolution[i].octave); 195 int sigma_size_ = fRound(evolution[i].esigma * options_.derivative_factor / ratio); 196 197 compute_scharr_derivatives(evolution[i].Lsmooth, evolution[i].Lx, 1, 0, sigma_size_); 198 compute_scharr_derivatives(evolution[i].Lsmooth, evolution[i].Ly, 0, 1, sigma_size_); 199 compute_scharr_derivatives(evolution[i].Lx, evolution[i].Lxx, 1, 0, sigma_size_); 200 compute_scharr_derivatives(evolution[i].Ly, evolution[i].Lyy, 0, 1, sigma_size_); 201 compute_scharr_derivatives(evolution[i].Lx, evolution[i].Lxy, 0, 1, sigma_size_); 202 203 evolution[i].Lx = evolution[i].Lx*((sigma_size_)); 204 evolution[i].Ly = evolution[i].Ly*((sigma_size_)); 205 evolution[i].Lxx = evolution[i].Lxx*((sigma_size_)*(sigma_size_)); 206 evolution[i].Lxy = evolution[i].Lxy*((sigma_size_)*(sigma_size_)); 207 evolution[i].Lyy = evolution[i].Lyy*((sigma_size_)*(sigma_size_)); 208 } 209 } 210 211 private: 212 std::vector<TEvolution>* evolution_; 213 AKAZEOptions options_; 214 }; 215 216 /* ************************************************************************* */ 217 /** 218 * @brief This method computes the multiscale derivatives for the nonlinear scale space 219 */ 220 void AKAZEFeatures::Compute_Multiscale_Derivatives(void) 221 { 222 parallel_for_(Range(0, (int)evolution_.size()), 223 MultiscaleDerivativesAKAZEInvoker(evolution_, options_)); 224 } 225 226 /* ************************************************************************* */ 227 /** 228 * @brief This method computes the feature detector response for the nonlinear scale space 229 * @note We use the Hessian determinant as the feature detector response 230 */ 231 void AKAZEFeatures::Compute_Determinant_Hessian_Response(void) { 232 233 // Firstly compute the multiscale derivatives 234 Compute_Multiscale_Derivatives(); 235 236 for (size_t i = 0; i < evolution_.size(); i++) 237 { 238 for (int ix = 0; ix < evolution_[i].Ldet.rows; ix++) 239 { 240 for (int jx = 0; jx < evolution_[i].Ldet.cols; jx++) 241 { 242 float lxx = *(evolution_[i].Lxx.ptr<float>(ix)+jx); 243 float lxy = *(evolution_[i].Lxy.ptr<float>(ix)+jx); 244 float lyy = *(evolution_[i].Lyy.ptr<float>(ix)+jx); 245 *(evolution_[i].Ldet.ptr<float>(ix)+jx) = (lxx*lyy - lxy*lxy); 246 } 247 } 248 } 249 } 250 251 /* ************************************************************************* */ 252 /** 253 * @brief This method finds extrema in the nonlinear scale space 254 * @param kpts Vector of detected keypoints 255 */ 256 void AKAZEFeatures::Find_Scale_Space_Extrema(std::vector<KeyPoint>& kpts) 257 { 258 259 float value = 0.0; 260 float dist = 0.0, ratio = 0.0, smax = 0.0; 261 int npoints = 0, id_repeated = 0; 262 int sigma_size_ = 0, left_x = 0, right_x = 0, up_y = 0, down_y = 0; 263 bool is_extremum = false, is_repeated = false, is_out = false; 264 KeyPoint point; 265 vector<KeyPoint> kpts_aux; 266 267 // Set maximum size 268 if (options_.descriptor == AKAZE::DESCRIPTOR_MLDB_UPRIGHT || options_.descriptor == AKAZE::DESCRIPTOR_MLDB) { 269 smax = 10.0f*sqrtf(2.0f); 270 } 271 else if (options_.descriptor == AKAZE::DESCRIPTOR_KAZE_UPRIGHT || options_.descriptor == AKAZE::DESCRIPTOR_KAZE) { 272 smax = 12.0f*sqrtf(2.0f); 273 } 274 275 for (size_t i = 0; i < evolution_.size(); i++) { 276 float* prev = evolution_[i].Ldet.ptr<float>(0); 277 float* curr = evolution_[i].Ldet.ptr<float>(1); 278 for (int ix = 1; ix < evolution_[i].Ldet.rows - 1; ix++) { 279 float* next = evolution_[i].Ldet.ptr<float>(ix + 1); 280 281 for (int jx = 1; jx < evolution_[i].Ldet.cols - 1; jx++) { 282 is_extremum = false; 283 is_repeated = false; 284 is_out = false; 285 value = *(evolution_[i].Ldet.ptr<float>(ix)+jx); 286 287 // Filter the points with the detector threshold 288 if (value > options_.dthreshold && value >= options_.min_dthreshold && 289 value > curr[jx-1] && 290 value > curr[jx+1] && 291 value > prev[jx-1] && 292 value > prev[jx] && 293 value > prev[jx+1] && 294 value > next[jx-1] && 295 value > next[jx] && 296 value > next[jx+1]) { 297 298 is_extremum = true; 299 point.response = fabs(value); 300 point.size = evolution_[i].esigma*options_.derivative_factor; 301 point.octave = (int)evolution_[i].octave; 302 point.class_id = (int)i; 303 ratio = (float)fastpow(2, point.octave); 304 sigma_size_ = fRound(point.size / ratio); 305 point.pt.x = static_cast<float>(jx); 306 point.pt.y = static_cast<float>(ix); 307 308 // Compare response with the same and lower scale 309 for (size_t ik = 0; ik < kpts_aux.size(); ik++) { 310 311 if ((point.class_id - 1) == kpts_aux[ik].class_id || 312 point.class_id == kpts_aux[ik].class_id) { 313 float distx = point.pt.x*ratio - kpts_aux[ik].pt.x; 314 float disty = point.pt.y*ratio - kpts_aux[ik].pt.y; 315 dist = distx * distx + disty * disty; 316 if (dist <= point.size * point.size) { 317 if (point.response > kpts_aux[ik].response) { 318 id_repeated = (int)ik; 319 is_repeated = true; 320 } 321 else { 322 is_extremum = false; 323 } 324 break; 325 } 326 } 327 } 328 329 // Check out of bounds 330 if (is_extremum == true) { 331 332 // Check that the point is under the image limits for the descriptor computation 333 left_x = fRound(point.pt.x - smax*sigma_size_) - 1; 334 right_x = fRound(point.pt.x + smax*sigma_size_) + 1; 335 up_y = fRound(point.pt.y - smax*sigma_size_) - 1; 336 down_y = fRound(point.pt.y + smax*sigma_size_) + 1; 337 338 if (left_x < 0 || right_x >= evolution_[i].Ldet.cols || 339 up_y < 0 || down_y >= evolution_[i].Ldet.rows) { 340 is_out = true; 341 } 342 343 if (is_out == false) { 344 if (is_repeated == false) { 345 point.pt.x *= ratio; 346 point.pt.y *= ratio; 347 kpts_aux.push_back(point); 348 npoints++; 349 } 350 else { 351 point.pt.x *= ratio; 352 point.pt.y *= ratio; 353 kpts_aux[id_repeated] = point; 354 } 355 } // if is_out 356 } //if is_extremum 357 } 358 } // for jx 359 prev = curr; 360 curr = next; 361 } // for ix 362 } // for i 363 364 // Now filter points with the upper scale level 365 for (size_t i = 0; i < kpts_aux.size(); i++) { 366 367 is_repeated = false; 368 const KeyPoint& pt = kpts_aux[i]; 369 for (size_t j = i + 1; j < kpts_aux.size(); j++) { 370 371 // Compare response with the upper scale 372 if ((pt.class_id + 1) == kpts_aux[j].class_id) { 373 float distx = pt.pt.x - kpts_aux[j].pt.x; 374 float disty = pt.pt.y - kpts_aux[j].pt.y; 375 dist = distx * distx + disty * disty; 376 if (dist <= pt.size * pt.size) { 377 if (pt.response < kpts_aux[j].response) { 378 is_repeated = true; 379 break; 380 } 381 } 382 } 383 } 384 385 if (is_repeated == false) 386 kpts.push_back(pt); 387 } 388 } 389 390 /* ************************************************************************* */ 391 /** 392 * @brief This method performs subpixel refinement of the detected keypoints 393 * @param kpts Vector of detected keypoints 394 */ 395 void AKAZEFeatures::Do_Subpixel_Refinement(std::vector<KeyPoint>& kpts) 396 { 397 float Dx = 0.0, Dy = 0.0, ratio = 0.0; 398 float Dxx = 0.0, Dyy = 0.0, Dxy = 0.0; 399 int x = 0, y = 0; 400 Matx22f A(0, 0, 0, 0); 401 Vec2f b(0, 0); 402 Vec2f dst(0, 0); 403 404 for (size_t i = 0; i < kpts.size(); i++) { 405 ratio = (float)fastpow(2, kpts[i].octave); 406 x = fRound(kpts[i].pt.x / ratio); 407 y = fRound(kpts[i].pt.y / ratio); 408 409 // Compute the gradient 410 Dx = (0.5f)*(*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y)+x + 1) 411 - *(evolution_[kpts[i].class_id].Ldet.ptr<float>(y)+x - 1)); 412 Dy = (0.5f)*(*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y + 1) + x) 413 - *(evolution_[kpts[i].class_id].Ldet.ptr<float>(y - 1) + x)); 414 415 // Compute the Hessian 416 Dxx = (*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y)+x + 1) 417 + *(evolution_[kpts[i].class_id].Ldet.ptr<float>(y)+x - 1) 418 - 2.0f*(*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y)+x))); 419 420 Dyy = (*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y + 1) + x) 421 + *(evolution_[kpts[i].class_id].Ldet.ptr<float>(y - 1) + x) 422 - 2.0f*(*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y)+x))); 423 424 Dxy = (0.25f)*(*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y + 1) + x + 1) 425 + (*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y - 1) + x - 1))) 426 - (0.25f)*(*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y - 1) + x + 1) 427 + (*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y + 1) + x - 1))); 428 429 // Solve the linear system 430 A(0, 0) = Dxx; 431 A(1, 1) = Dyy; 432 A(0, 1) = A(1, 0) = Dxy; 433 b(0) = -Dx; 434 b(1) = -Dy; 435 436 solve(A, b, dst, DECOMP_LU); 437 438 if (fabs(dst(0)) <= 1.0f && fabs(dst(1)) <= 1.0f) { 439 kpts[i].pt.x = x + dst(0); 440 kpts[i].pt.y = y + dst(1); 441 int power = fastpow(2, evolution_[kpts[i].class_id].octave); 442 kpts[i].pt.x *= power; 443 kpts[i].pt.y *= power; 444 kpts[i].angle = 0.0; 445 446 // In OpenCV the size of a keypoint its the diameter 447 kpts[i].size *= 2.0f; 448 } 449 // Delete the point since its not stable 450 else { 451 kpts.erase(kpts.begin() + i); 452 i--; 453 } 454 } 455 } 456 457 /* ************************************************************************* */ 458 459 class SURF_Descriptor_Upright_64_Invoker : public ParallelLoopBody 460 { 461 public: 462 SURF_Descriptor_Upright_64_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, std::vector<TEvolution>& evolution) 463 : keypoints_(&kpts) 464 , descriptors_(&desc) 465 , evolution_(&evolution) 466 { 467 } 468 469 void operator() (const Range& range) const 470 { 471 for (int i = range.start; i < range.end; i++) 472 { 473 Get_SURF_Descriptor_Upright_64((*keypoints_)[i], descriptors_->ptr<float>(i)); 474 } 475 } 476 477 void Get_SURF_Descriptor_Upright_64(const KeyPoint& kpt, float* desc) const; 478 479 private: 480 std::vector<KeyPoint>* keypoints_; 481 Mat* descriptors_; 482 std::vector<TEvolution>* evolution_; 483 }; 484 485 class SURF_Descriptor_64_Invoker : public ParallelLoopBody 486 { 487 public: 488 SURF_Descriptor_64_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, std::vector<TEvolution>& evolution) 489 : keypoints_(&kpts) 490 , descriptors_(&desc) 491 , evolution_(&evolution) 492 { 493 } 494 495 void operator()(const Range& range) const 496 { 497 for (int i = range.start; i < range.end; i++) 498 { 499 AKAZEFeatures::Compute_Main_Orientation((*keypoints_)[i], *evolution_); 500 Get_SURF_Descriptor_64((*keypoints_)[i], descriptors_->ptr<float>(i)); 501 } 502 } 503 504 void Get_SURF_Descriptor_64(const KeyPoint& kpt, float* desc) const; 505 506 private: 507 std::vector<KeyPoint>* keypoints_; 508 Mat* descriptors_; 509 std::vector<TEvolution>* evolution_; 510 }; 511 512 class MSURF_Upright_Descriptor_64_Invoker : public ParallelLoopBody 513 { 514 public: 515 MSURF_Upright_Descriptor_64_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, std::vector<TEvolution>& evolution) 516 : keypoints_(&kpts) 517 , descriptors_(&desc) 518 , evolution_(&evolution) 519 { 520 } 521 522 void operator()(const Range& range) const 523 { 524 for (int i = range.start; i < range.end; i++) 525 { 526 Get_MSURF_Upright_Descriptor_64((*keypoints_)[i], descriptors_->ptr<float>(i)); 527 } 528 } 529 530 void Get_MSURF_Upright_Descriptor_64(const KeyPoint& kpt, float* desc) const; 531 532 private: 533 std::vector<KeyPoint>* keypoints_; 534 Mat* descriptors_; 535 std::vector<TEvolution>* evolution_; 536 }; 537 538 class MSURF_Descriptor_64_Invoker : public ParallelLoopBody 539 { 540 public: 541 MSURF_Descriptor_64_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, std::vector<TEvolution>& evolution) 542 : keypoints_(&kpts) 543 , descriptors_(&desc) 544 , evolution_(&evolution) 545 { 546 } 547 548 void operator() (const Range& range) const 549 { 550 for (int i = range.start; i < range.end; i++) 551 { 552 AKAZEFeatures::Compute_Main_Orientation((*keypoints_)[i], *evolution_); 553 Get_MSURF_Descriptor_64((*keypoints_)[i], descriptors_->ptr<float>(i)); 554 } 555 } 556 557 void Get_MSURF_Descriptor_64(const KeyPoint& kpt, float* desc) const; 558 559 private: 560 std::vector<KeyPoint>* keypoints_; 561 Mat* descriptors_; 562 std::vector<TEvolution>* evolution_; 563 }; 564 565 class Upright_MLDB_Full_Descriptor_Invoker : public ParallelLoopBody 566 { 567 public: 568 Upright_MLDB_Full_Descriptor_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, std::vector<TEvolution>& evolution, AKAZEOptions& options) 569 : keypoints_(&kpts) 570 , descriptors_(&desc) 571 , evolution_(&evolution) 572 , options_(&options) 573 { 574 } 575 576 void operator() (const Range& range) const 577 { 578 for (int i = range.start; i < range.end; i++) 579 { 580 Get_Upright_MLDB_Full_Descriptor((*keypoints_)[i], descriptors_->ptr<unsigned char>(i)); 581 } 582 } 583 584 void Get_Upright_MLDB_Full_Descriptor(const KeyPoint& kpt, unsigned char* desc) const; 585 586 private: 587 std::vector<KeyPoint>* keypoints_; 588 Mat* descriptors_; 589 std::vector<TEvolution>* evolution_; 590 AKAZEOptions* options_; 591 }; 592 593 class Upright_MLDB_Descriptor_Subset_Invoker : public ParallelLoopBody 594 { 595 public: 596 Upright_MLDB_Descriptor_Subset_Invoker(std::vector<KeyPoint>& kpts, 597 Mat& desc, 598 std::vector<TEvolution>& evolution, 599 AKAZEOptions& options, 600 Mat descriptorSamples, 601 Mat descriptorBits) 602 : keypoints_(&kpts) 603 , descriptors_(&desc) 604 , evolution_(&evolution) 605 , options_(&options) 606 , descriptorSamples_(descriptorSamples) 607 , descriptorBits_(descriptorBits) 608 { 609 } 610 611 void operator() (const Range& range) const 612 { 613 for (int i = range.start; i < range.end; i++) 614 { 615 Get_Upright_MLDB_Descriptor_Subset((*keypoints_)[i], descriptors_->ptr<unsigned char>(i)); 616 } 617 } 618 619 void Get_Upright_MLDB_Descriptor_Subset(const KeyPoint& kpt, unsigned char* desc) const; 620 621 private: 622 std::vector<KeyPoint>* keypoints_; 623 Mat* descriptors_; 624 std::vector<TEvolution>* evolution_; 625 AKAZEOptions* options_; 626 627 Mat descriptorSamples_; // List of positions in the grids to sample LDB bits from. 628 Mat descriptorBits_; 629 }; 630 631 class MLDB_Full_Descriptor_Invoker : public ParallelLoopBody 632 { 633 public: 634 MLDB_Full_Descriptor_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, std::vector<TEvolution>& evolution, AKAZEOptions& options) 635 : keypoints_(&kpts) 636 , descriptors_(&desc) 637 , evolution_(&evolution) 638 , options_(&options) 639 { 640 } 641 642 void operator() (const Range& range) const 643 { 644 for (int i = range.start; i < range.end; i++) 645 { 646 AKAZEFeatures::Compute_Main_Orientation((*keypoints_)[i], *evolution_); 647 Get_MLDB_Full_Descriptor((*keypoints_)[i], descriptors_->ptr<unsigned char>(i)); 648 } 649 } 650 651 void Get_MLDB_Full_Descriptor(const KeyPoint& kpt, unsigned char* desc) const; 652 void MLDB_Fill_Values(float* values, int sample_step, int level, 653 float xf, float yf, float co, float si, float scale) const; 654 void MLDB_Binary_Comparisons(float* values, unsigned char* desc, 655 int count, int& dpos) const; 656 657 private: 658 std::vector<KeyPoint>* keypoints_; 659 Mat* descriptors_; 660 std::vector<TEvolution>* evolution_; 661 AKAZEOptions* options_; 662 }; 663 664 class MLDB_Descriptor_Subset_Invoker : public ParallelLoopBody 665 { 666 public: 667 MLDB_Descriptor_Subset_Invoker(std::vector<KeyPoint>& kpts, 668 Mat& desc, 669 std::vector<TEvolution>& evolution, 670 AKAZEOptions& options, 671 Mat descriptorSamples, 672 Mat descriptorBits) 673 : keypoints_(&kpts) 674 , descriptors_(&desc) 675 , evolution_(&evolution) 676 , options_(&options) 677 , descriptorSamples_(descriptorSamples) 678 , descriptorBits_(descriptorBits) 679 { 680 } 681 682 void operator() (const Range& range) const 683 { 684 for (int i = range.start; i < range.end; i++) 685 { 686 AKAZEFeatures::Compute_Main_Orientation((*keypoints_)[i], *evolution_); 687 Get_MLDB_Descriptor_Subset((*keypoints_)[i], descriptors_->ptr<unsigned char>(i)); 688 } 689 } 690 691 void Get_MLDB_Descriptor_Subset(const KeyPoint& kpt, unsigned char* desc) const; 692 693 private: 694 std::vector<KeyPoint>* keypoints_; 695 Mat* descriptors_; 696 std::vector<TEvolution>* evolution_; 697 AKAZEOptions* options_; 698 699 Mat descriptorSamples_; // List of positions in the grids to sample LDB bits from. 700 Mat descriptorBits_; 701 }; 702 703 /** 704 * @brief This method computes the set of descriptors through the nonlinear scale space 705 * @param kpts Vector of detected keypoints 706 * @param desc Matrix to store the descriptors 707 */ 708 void AKAZEFeatures::Compute_Descriptors(std::vector<KeyPoint>& kpts, Mat& desc) 709 { 710 for(size_t i = 0; i < kpts.size(); i++) 711 { 712 CV_Assert(0 <= kpts[i].class_id && kpts[i].class_id < static_cast<int>(evolution_.size())); 713 } 714 715 // Allocate memory for the matrix with the descriptors 716 if (options_.descriptor < AKAZE::DESCRIPTOR_MLDB_UPRIGHT) { 717 desc = Mat::zeros((int)kpts.size(), 64, CV_32FC1); 718 } 719 else { 720 // We use the full length binary descriptor -> 486 bits 721 if (options_.descriptor_size == 0) { 722 int t = (6 + 36 + 120)*options_.descriptor_channels; 723 desc = Mat::zeros((int)kpts.size(), (int)ceil(t / 8.), CV_8UC1); 724 } 725 else { 726 // We use the random bit selection length binary descriptor 727 desc = Mat::zeros((int)kpts.size(), (int)ceil(options_.descriptor_size / 8.), CV_8UC1); 728 } 729 } 730 731 switch (options_.descriptor) 732 { 733 case AKAZE::DESCRIPTOR_KAZE_UPRIGHT: // Upright descriptors, not invariant to rotation 734 { 735 parallel_for_(Range(0, (int)kpts.size()), MSURF_Upright_Descriptor_64_Invoker(kpts, desc, evolution_)); 736 } 737 break; 738 case AKAZE::DESCRIPTOR_KAZE: 739 { 740 parallel_for_(Range(0, (int)kpts.size()), MSURF_Descriptor_64_Invoker(kpts, desc, evolution_)); 741 } 742 break; 743 case AKAZE::DESCRIPTOR_MLDB_UPRIGHT: // Upright descriptors, not invariant to rotation 744 { 745 if (options_.descriptor_size == 0) 746 parallel_for_(Range(0, (int)kpts.size()), Upright_MLDB_Full_Descriptor_Invoker(kpts, desc, evolution_, options_)); 747 else 748 parallel_for_(Range(0, (int)kpts.size()), Upright_MLDB_Descriptor_Subset_Invoker(kpts, desc, evolution_, options_, descriptorSamples_, descriptorBits_)); 749 } 750 break; 751 case AKAZE::DESCRIPTOR_MLDB: 752 { 753 if (options_.descriptor_size == 0) 754 parallel_for_(Range(0, (int)kpts.size()), MLDB_Full_Descriptor_Invoker(kpts, desc, evolution_, options_)); 755 else 756 parallel_for_(Range(0, (int)kpts.size()), MLDB_Descriptor_Subset_Invoker(kpts, desc, evolution_, options_, descriptorSamples_, descriptorBits_)); 757 } 758 break; 759 } 760 } 761 762 /* ************************************************************************* */ 763 /** 764 * @brief This method computes the main orientation for a given keypoint 765 * @param kpt Input keypoint 766 * @note The orientation is computed using a similar approach as described in the 767 * original SURF method. See Bay et al., Speeded Up Robust Features, ECCV 2006 768 */ 769 void AKAZEFeatures::Compute_Main_Orientation(KeyPoint& kpt, const std::vector<TEvolution>& evolution_) 770 { 771 /* ************************************************************************* */ 772 /// Lookup table for 2d gaussian (sigma = 2.5) where (0,0) is top left and (6,6) is bottom right 773 static const float gauss25[7][7] = 774 { 775 { 0.02546481f, 0.02350698f, 0.01849125f, 0.01239505f, 0.00708017f, 0.00344629f, 0.00142946f }, 776 { 0.02350698f, 0.02169968f, 0.01706957f, 0.01144208f, 0.00653582f, 0.00318132f, 0.00131956f }, 777 { 0.01849125f, 0.01706957f, 0.01342740f, 0.00900066f, 0.00514126f, 0.00250252f, 0.00103800f }, 778 { 0.01239505f, 0.01144208f, 0.00900066f, 0.00603332f, 0.00344629f, 0.00167749f, 0.00069579f }, 779 { 0.00708017f, 0.00653582f, 0.00514126f, 0.00344629f, 0.00196855f, 0.00095820f, 0.00039744f }, 780 { 0.00344629f, 0.00318132f, 0.00250252f, 0.00167749f, 0.00095820f, 0.00046640f, 0.00019346f }, 781 { 0.00142946f, 0.00131956f, 0.00103800f, 0.00069579f, 0.00039744f, 0.00019346f, 0.00008024f } 782 }; 783 784 int ix = 0, iy = 0, idx = 0, s = 0, level = 0; 785 float xf = 0.0, yf = 0.0, gweight = 0.0, ratio = 0.0; 786 const int ang_size = 109; 787 float resX[ang_size], resY[ang_size], Ang[ang_size]; 788 const int id[] = { 6, 5, 4, 3, 2, 1, 0, 1, 2, 3, 4, 5, 6 }; 789 790 // Variables for computing the dominant direction 791 float sumX = 0.0, sumY = 0.0, max = 0.0, ang1 = 0.0, ang2 = 0.0; 792 793 // Get the information from the keypoint 794 level = kpt.class_id; 795 ratio = (float)(1 << evolution_[level].octave); 796 s = fRound(0.5f*kpt.size / ratio); 797 xf = kpt.pt.x / ratio; 798 yf = kpt.pt.y / ratio; 799 800 // Calculate derivatives responses for points within radius of 6*scale 801 for (int i = -6; i <= 6; ++i) { 802 for (int j = -6; j <= 6; ++j) { 803 if (i*i + j*j < 36) { 804 iy = fRound(yf + j*s); 805 ix = fRound(xf + i*s); 806 807 gweight = gauss25[id[i + 6]][id[j + 6]]; 808 resX[idx] = gweight*(*(evolution_[level].Lx.ptr<float>(iy)+ix)); 809 resY[idx] = gweight*(*(evolution_[level].Ly.ptr<float>(iy)+ix)); 810 811 ++idx; 812 } 813 } 814 } 815 hal::fastAtan2(resY, resX, Ang, ang_size, false); 816 // Loop slides pi/3 window around feature point 817 for (ang1 = 0; ang1 < (float)(2.0 * CV_PI); ang1 += 0.15f) { 818 ang2 = (ang1 + (float)(CV_PI / 3.0) >(float)(2.0*CV_PI) ? ang1 - (float)(5.0*CV_PI / 3.0) : ang1 + (float)(CV_PI / 3.0)); 819 sumX = sumY = 0.f; 820 821 for (int k = 0; k < ang_size; ++k) { 822 // Get angle from the x-axis of the sample point 823 const float & ang = Ang[k]; 824 825 // Determine whether the point is within the window 826 if (ang1 < ang2 && ang1 < ang && ang < ang2) { 827 sumX += resX[k]; 828 sumY += resY[k]; 829 } 830 else if (ang2 < ang1 && 831 ((ang > 0 && ang < ang2) || (ang > ang1 && ang < 2.0f*CV_PI))) { 832 sumX += resX[k]; 833 sumY += resY[k]; 834 } 835 } 836 837 // if the vector produced from this window is longer than all 838 // previous vectors then this forms the new dominant direction 839 if (sumX*sumX + sumY*sumY > max) { 840 // store largest orientation 841 max = sumX*sumX + sumY*sumY; 842 kpt.angle = getAngle(sumX, sumY); 843 } 844 } 845 } 846 847 /* ************************************************************************* */ 848 /** 849 * @brief This method computes the upright descriptor (not rotation invariant) of 850 * the provided keypoint 851 * @param kpt Input keypoint 852 * @param desc Descriptor vector 853 * @note Rectangular grid of 24 s x 24 s. Descriptor Length 64. The descriptor is inspired 854 * from Agrawal et al., CenSurE: Center Surround Extremas for Realtime Feature Detection and Matching, 855 * ECCV 2008 856 */ 857 void MSURF_Upright_Descriptor_64_Invoker::Get_MSURF_Upright_Descriptor_64(const KeyPoint& kpt, float *desc) const { 858 859 float dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0, gauss_s1 = 0.0, gauss_s2 = 0.0; 860 float rx = 0.0, ry = 0.0, len = 0.0, xf = 0.0, yf = 0.0, ys = 0.0, xs = 0.0; 861 float sample_x = 0.0, sample_y = 0.0; 862 int x1 = 0, y1 = 0, sample_step = 0, pattern_size = 0; 863 int x2 = 0, y2 = 0, kx = 0, ky = 0, i = 0, j = 0, dcount = 0; 864 float fx = 0.0, fy = 0.0, ratio = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0, res4 = 0.0; 865 int scale = 0, dsize = 0, level = 0; 866 867 // Subregion centers for the 4x4 gaussian weighting 868 float cx = -0.5f, cy = 0.5f; 869 870 const std::vector<TEvolution>& evolution = *evolution_; 871 872 // Set the descriptor size and the sample and pattern sizes 873 dsize = 64; 874 sample_step = 5; 875 pattern_size = 12; 876 877 // Get the information from the keypoint 878 ratio = (float)(1 << kpt.octave); 879 scale = fRound(0.5f*kpt.size / ratio); 880 level = kpt.class_id; 881 yf = kpt.pt.y / ratio; 882 xf = kpt.pt.x / ratio; 883 884 i = -8; 885 886 // Calculate descriptor for this interest point 887 // Area of size 24 s x 24 s 888 while (i < pattern_size) { 889 j = -8; 890 i = i - 4; 891 892 cx += 1.0f; 893 cy = -0.5f; 894 895 while (j < pattern_size) { 896 dx = dy = mdx = mdy = 0.0; 897 cy += 1.0f; 898 j = j - 4; 899 900 ky = i + sample_step; 901 kx = j + sample_step; 902 903 ys = yf + (ky*scale); 904 xs = xf + (kx*scale); 905 906 for (int k = i; k < i + 9; k++) { 907 for (int l = j; l < j + 9; l++) { 908 sample_y = k*scale + yf; 909 sample_x = l*scale + xf; 910 911 //Get the gaussian weighted x and y responses 912 gauss_s1 = gaussian(xs - sample_x, ys - sample_y, 2.50f*scale); 913 914 y1 = (int)(sample_y - .5); 915 x1 = (int)(sample_x - .5); 916 917 y2 = (int)(sample_y + .5); 918 x2 = (int)(sample_x + .5); 919 920 fx = sample_x - x1; 921 fy = sample_y - y1; 922 923 res1 = *(evolution[level].Lx.ptr<float>(y1)+x1); 924 res2 = *(evolution[level].Lx.ptr<float>(y1)+x2); 925 res3 = *(evolution[level].Lx.ptr<float>(y2)+x1); 926 res4 = *(evolution[level].Lx.ptr<float>(y2)+x2); 927 rx = (1.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + (1.0f - fx)*fy*res3 + fx*fy*res4; 928 929 res1 = *(evolution[level].Ly.ptr<float>(y1)+x1); 930 res2 = *(evolution[level].Ly.ptr<float>(y1)+x2); 931 res3 = *(evolution[level].Ly.ptr<float>(y2)+x1); 932 res4 = *(evolution[level].Ly.ptr<float>(y2)+x2); 933 ry = (1.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + (1.0f - fx)*fy*res3 + fx*fy*res4; 934 935 rx = gauss_s1*rx; 936 ry = gauss_s1*ry; 937 938 // Sum the derivatives to the cumulative descriptor 939 dx += rx; 940 dy += ry; 941 mdx += fabs(rx); 942 mdy += fabs(ry); 943 } 944 } 945 946 // Add the values to the descriptor vector 947 gauss_s2 = gaussian(cx - 2.0f, cy - 2.0f, 1.5f); 948 949 desc[dcount++] = dx*gauss_s2; 950 desc[dcount++] = dy*gauss_s2; 951 desc[dcount++] = mdx*gauss_s2; 952 desc[dcount++] = mdy*gauss_s2; 953 954 len += (dx*dx + dy*dy + mdx*mdx + mdy*mdy)*gauss_s2*gauss_s2; 955 956 j += 9; 957 } 958 959 i += 9; 960 } 961 962 // convert to unit vector 963 len = sqrt(len); 964 965 for (i = 0; i < dsize; i++) { 966 desc[i] /= len; 967 } 968 } 969 970 /* ************************************************************************* */ 971 /** 972 * @brief This method computes the descriptor of the provided keypoint given the 973 * main orientation of the keypoint 974 * @param kpt Input keypoint 975 * @param desc Descriptor vector 976 * @note Rectangular grid of 24 s x 24 s. Descriptor Length 64. The descriptor is inspired 977 * from Agrawal et al., CenSurE: Center Surround Extremas for Realtime Feature Detection and Matching, 978 * ECCV 2008 979 */ 980 void MSURF_Descriptor_64_Invoker::Get_MSURF_Descriptor_64(const KeyPoint& kpt, float *desc) const { 981 982 float dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0, gauss_s1 = 0.0, gauss_s2 = 0.0; 983 float rx = 0.0, ry = 0.0, rrx = 0.0, rry = 0.0, len = 0.0, xf = 0.0, yf = 0.0, ys = 0.0, xs = 0.0; 984 float sample_x = 0.0, sample_y = 0.0, co = 0.0, si = 0.0, angle = 0.0; 985 float fx = 0.0, fy = 0.0, ratio = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0, res4 = 0.0; 986 int x1 = 0, y1 = 0, x2 = 0, y2 = 0, sample_step = 0, pattern_size = 0; 987 int kx = 0, ky = 0, i = 0, j = 0, dcount = 0; 988 int scale = 0, dsize = 0, level = 0; 989 990 // Subregion centers for the 4x4 gaussian weighting 991 float cx = -0.5f, cy = 0.5f; 992 993 const std::vector<TEvolution>& evolution = *evolution_; 994 995 // Set the descriptor size and the sample and pattern sizes 996 dsize = 64; 997 sample_step = 5; 998 pattern_size = 12; 999 1000 // Get the information from the keypoint 1001 ratio = (float)(1 << kpt.octave); 1002 scale = fRound(0.5f*kpt.size / ratio); 1003 angle = kpt.angle; 1004 level = kpt.class_id; 1005 yf = kpt.pt.y / ratio; 1006 xf = kpt.pt.x / ratio; 1007 co = cos(angle); 1008 si = sin(angle); 1009 1010 i = -8; 1011 1012 // Calculate descriptor for this interest point 1013 // Area of size 24 s x 24 s 1014 while (i < pattern_size) { 1015 j = -8; 1016 i = i - 4; 1017 1018 cx += 1.0f; 1019 cy = -0.5f; 1020 1021 while (j < pattern_size) { 1022 dx = dy = mdx = mdy = 0.0; 1023 cy += 1.0f; 1024 j = j - 4; 1025 1026 ky = i + sample_step; 1027 kx = j + sample_step; 1028 1029 xs = xf + (-kx*scale*si + ky*scale*co); 1030 ys = yf + (kx*scale*co + ky*scale*si); 1031 1032 for (int k = i; k < i + 9; ++k) { 1033 for (int l = j; l < j + 9; ++l) { 1034 // Get coords of sample point on the rotated axis 1035 sample_y = yf + (l*scale*co + k*scale*si); 1036 sample_x = xf + (-l*scale*si + k*scale*co); 1037 1038 // Get the gaussian weighted x and y responses 1039 gauss_s1 = gaussian(xs - sample_x, ys - sample_y, 2.5f*scale); 1040 1041 y1 = fRound(sample_y - 0.5f); 1042 x1 = fRound(sample_x - 0.5f); 1043 1044 y2 = fRound(sample_y + 0.5f); 1045 x2 = fRound(sample_x + 0.5f); 1046 1047 fx = sample_x - x1; 1048 fy = sample_y - y1; 1049 1050 res1 = *(evolution[level].Lx.ptr<float>(y1)+x1); 1051 res2 = *(evolution[level].Lx.ptr<float>(y1)+x2); 1052 res3 = *(evolution[level].Lx.ptr<float>(y2)+x1); 1053 res4 = *(evolution[level].Lx.ptr<float>(y2)+x2); 1054 rx = (1.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + (1.0f - fx)*fy*res3 + fx*fy*res4; 1055 1056 res1 = *(evolution[level].Ly.ptr<float>(y1)+x1); 1057 res2 = *(evolution[level].Ly.ptr<float>(y1)+x2); 1058 res3 = *(evolution[level].Ly.ptr<float>(y2)+x1); 1059 res4 = *(evolution[level].Ly.ptr<float>(y2)+x2); 1060 ry = (1.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + (1.0f - fx)*fy*res3 + fx*fy*res4; 1061 1062 // Get the x and y derivatives on the rotated axis 1063 rry = gauss_s1*(rx*co + ry*si); 1064 rrx = gauss_s1*(-rx*si + ry*co); 1065 1066 // Sum the derivatives to the cumulative descriptor 1067 dx += rrx; 1068 dy += rry; 1069 mdx += fabs(rrx); 1070 mdy += fabs(rry); 1071 } 1072 } 1073 1074 // Add the values to the descriptor vector 1075 gauss_s2 = gaussian(cx - 2.0f, cy - 2.0f, 1.5f); 1076 desc[dcount++] = dx*gauss_s2; 1077 desc[dcount++] = dy*gauss_s2; 1078 desc[dcount++] = mdx*gauss_s2; 1079 desc[dcount++] = mdy*gauss_s2; 1080 1081 len += (dx*dx + dy*dy + mdx*mdx + mdy*mdy)*gauss_s2*gauss_s2; 1082 1083 j += 9; 1084 } 1085 1086 i += 9; 1087 } 1088 1089 // convert to unit vector 1090 len = sqrt(len); 1091 1092 for (i = 0; i < dsize; i++) { 1093 desc[i] /= len; 1094 } 1095 } 1096 1097 /* ************************************************************************* */ 1098 /** 1099 * @brief This method computes the rupright descriptor (not rotation invariant) of 1100 * the provided keypoint 1101 * @param kpt Input keypoint 1102 * @param desc Descriptor vector 1103 */ 1104 void Upright_MLDB_Full_Descriptor_Invoker::Get_Upright_MLDB_Full_Descriptor(const KeyPoint& kpt, unsigned char *desc) const { 1105 1106 float di = 0.0, dx = 0.0, dy = 0.0; 1107 float ri = 0.0, rx = 0.0, ry = 0.0, xf = 0.0, yf = 0.0; 1108 float sample_x = 0.0, sample_y = 0.0, ratio = 0.0; 1109 int x1 = 0, y1 = 0, sample_step = 0, pattern_size = 0; 1110 int level = 0, nsamples = 0, scale = 0; 1111 int dcount1 = 0, dcount2 = 0; 1112 1113 const AKAZEOptions & options = *options_; 1114 const std::vector<TEvolution>& evolution = *evolution_; 1115 1116 // Matrices for the M-LDB descriptor 1117 Mat values_1 = Mat::zeros(4, options.descriptor_channels, CV_32FC1); 1118 Mat values_2 = Mat::zeros(9, options.descriptor_channels, CV_32FC1); 1119 Mat values_3 = Mat::zeros(16, options.descriptor_channels, CV_32FC1); 1120 1121 // Get the information from the keypoint 1122 ratio = (float)(1 << kpt.octave); 1123 scale = fRound(0.5f*kpt.size / ratio); 1124 level = kpt.class_id; 1125 yf = kpt.pt.y / ratio; 1126 xf = kpt.pt.x / ratio; 1127 1128 // First 2x2 grid 1129 pattern_size = options_->descriptor_pattern_size; 1130 sample_step = pattern_size; 1131 1132 for (int i = -pattern_size; i < pattern_size; i += sample_step) { 1133 for (int j = -pattern_size; j < pattern_size; j += sample_step) { 1134 di = dx = dy = 0.0; 1135 nsamples = 0; 1136 1137 for (int k = i; k < i + sample_step; k++) { 1138 for (int l = j; l < j + sample_step; l++) { 1139 1140 // Get the coordinates of the sample point 1141 sample_y = yf + l*scale; 1142 sample_x = xf + k*scale; 1143 1144 y1 = fRound(sample_y); 1145 x1 = fRound(sample_x); 1146 1147 ri = *(evolution[level].Lt.ptr<float>(y1)+x1); 1148 rx = *(evolution[level].Lx.ptr<float>(y1)+x1); 1149 ry = *(evolution[level].Ly.ptr<float>(y1)+x1); 1150 1151 di += ri; 1152 dx += rx; 1153 dy += ry; 1154 nsamples++; 1155 } 1156 } 1157 1158 di /= nsamples; 1159 dx /= nsamples; 1160 dy /= nsamples; 1161 1162 *(values_1.ptr<float>(dcount2)) = di; 1163 *(values_1.ptr<float>(dcount2)+1) = dx; 1164 *(values_1.ptr<float>(dcount2)+2) = dy; 1165 dcount2++; 1166 } 1167 } 1168 1169 // Do binary comparison first level 1170 for (int i = 0; i < 4; i++) { 1171 for (int j = i + 1; j < 4; j++) { 1172 if (*(values_1.ptr<float>(i)) > *(values_1.ptr<float>(j))) { 1173 desc[dcount1 / 8] |= (1 << (dcount1 % 8)); 1174 } 1175 dcount1++; 1176 1177 if (*(values_1.ptr<float>(i)+1) > *(values_1.ptr<float>(j)+1)) { 1178 desc[dcount1 / 8] |= (1 << (dcount1 % 8)); 1179 } 1180 dcount1++; 1181 1182 if (*(values_1.ptr<float>(i)+2) > *(values_1.ptr<float>(j)+2)) { 1183 desc[dcount1 / 8] |= (1 << (dcount1 % 8)); 1184 } 1185 dcount1++; 1186 } 1187 } 1188 1189 // Second 3x3 grid 1190 sample_step = static_cast<int>(ceil(pattern_size*2. / 3.)); 1191 dcount2 = 0; 1192 1193 for (int i = -pattern_size; i < pattern_size; i += sample_step) { 1194 for (int j = -pattern_size; j < pattern_size; j += sample_step) { 1195 di = dx = dy = 0.0; 1196 nsamples = 0; 1197 1198 for (int k = i; k < i + sample_step; k++) { 1199 for (int l = j; l < j + sample_step; l++) { 1200 1201 // Get the coordinates of the sample point 1202 sample_y = yf + l*scale; 1203 sample_x = xf + k*scale; 1204 1205 y1 = fRound(sample_y); 1206 x1 = fRound(sample_x); 1207 1208 ri = *(evolution[level].Lt.ptr<float>(y1)+x1); 1209 rx = *(evolution[level].Lx.ptr<float>(y1)+x1); 1210 ry = *(evolution[level].Ly.ptr<float>(y1)+x1); 1211 1212 di += ri; 1213 dx += rx; 1214 dy += ry; 1215 nsamples++; 1216 } 1217 } 1218 1219 di /= nsamples; 1220 dx /= nsamples; 1221 dy /= nsamples; 1222 1223 *(values_2.ptr<float>(dcount2)) = di; 1224 *(values_2.ptr<float>(dcount2)+1) = dx; 1225 *(values_2.ptr<float>(dcount2)+2) = dy; 1226 dcount2++; 1227 } 1228 } 1229 1230 //Do binary comparison second level 1231 dcount2 = 0; 1232 for (int i = 0; i < 9; i++) { 1233 for (int j = i + 1; j < 9; j++) { 1234 if (*(values_2.ptr<float>(i)) > *(values_2.ptr<float>(j))) { 1235 desc[dcount1 / 8] |= (1 << (dcount1 % 8)); 1236 } 1237 dcount1++; 1238 1239 if (*(values_2.ptr<float>(i)+1) > *(values_2.ptr<float>(j)+1)) { 1240 desc[dcount1 / 8] |= (1 << (dcount1 % 8)); 1241 } 1242 dcount1++; 1243 1244 if (*(values_2.ptr<float>(i)+2) > *(values_2.ptr<float>(j)+2)) { 1245 desc[dcount1 / 8] |= (1 << (dcount1 % 8)); 1246 } 1247 dcount1++; 1248 } 1249 } 1250 1251 // Third 4x4 grid 1252 sample_step = pattern_size / 2; 1253 dcount2 = 0; 1254 1255 for (int i = -pattern_size; i < pattern_size; i += sample_step) { 1256 for (int j = -pattern_size; j < pattern_size; j += sample_step) { 1257 di = dx = dy = 0.0; 1258 nsamples = 0; 1259 1260 for (int k = i; k < i + sample_step; k++) { 1261 for (int l = j; l < j + sample_step; l++) { 1262 1263 // Get the coordinates of the sample point 1264 sample_y = yf + l*scale; 1265 sample_x = xf + k*scale; 1266 1267 y1 = fRound(sample_y); 1268 x1 = fRound(sample_x); 1269 1270 ri = *(evolution[level].Lt.ptr<float>(y1)+x1); 1271 rx = *(evolution[level].Lx.ptr<float>(y1)+x1); 1272 ry = *(evolution[level].Ly.ptr<float>(y1)+x1); 1273 1274 di += ri; 1275 dx += rx; 1276 dy += ry; 1277 nsamples++; 1278 } 1279 } 1280 1281 di /= nsamples; 1282 dx /= nsamples; 1283 dy /= nsamples; 1284 1285 *(values_3.ptr<float>(dcount2)) = di; 1286 *(values_3.ptr<float>(dcount2)+1) = dx; 1287 *(values_3.ptr<float>(dcount2)+2) = dy; 1288 dcount2++; 1289 } 1290 } 1291 1292 //Do binary comparison third level 1293 dcount2 = 0; 1294 for (int i = 0; i < 16; i++) { 1295 for (int j = i + 1; j < 16; j++) { 1296 if (*(values_3.ptr<float>(i)) > *(values_3.ptr<float>(j))) { 1297 desc[dcount1 / 8] |= (1 << (dcount1 % 8)); 1298 } 1299 dcount1++; 1300 1301 if (*(values_3.ptr<float>(i)+1) > *(values_3.ptr<float>(j)+1)) { 1302 desc[dcount1 / 8] |= (1 << (dcount1 % 8)); 1303 } 1304 dcount1++; 1305 1306 if (*(values_3.ptr<float>(i)+2) > *(values_3.ptr<float>(j)+2)) { 1307 desc[dcount1 / 8] |= (1 << (dcount1 % 8)); 1308 } 1309 dcount1++; 1310 } 1311 } 1312 } 1313 1314 void MLDB_Full_Descriptor_Invoker::MLDB_Fill_Values(float* values, int sample_step, int level, 1315 float xf, float yf, float co, float si, float scale) const 1316 { 1317 const std::vector<TEvolution>& evolution = *evolution_; 1318 int pattern_size = options_->descriptor_pattern_size; 1319 int chan = options_->descriptor_channels; 1320 int valpos = 0; 1321 1322 for (int i = -pattern_size; i < pattern_size; i += sample_step) { 1323 for (int j = -pattern_size; j < pattern_size; j += sample_step) { 1324 float di, dx, dy; 1325 di = dx = dy = 0.0; 1326 int nsamples = 0; 1327 1328 for (int k = i; k < i + sample_step; k++) { 1329 for (int l = j; l < j + sample_step; l++) { 1330 float sample_y = yf + (l*co * scale + k*si*scale); 1331 float sample_x = xf + (-l*si * scale + k*co*scale); 1332 1333 int y1 = fRound(sample_y); 1334 int x1 = fRound(sample_x); 1335 1336 float ri = *(evolution[level].Lt.ptr<float>(y1)+x1); 1337 di += ri; 1338 1339 if(chan > 1) { 1340 float rx = *(evolution[level].Lx.ptr<float>(y1)+x1); 1341 float ry = *(evolution[level].Ly.ptr<float>(y1)+x1); 1342 if (chan == 2) { 1343 dx += sqrtf(rx*rx + ry*ry); 1344 } 1345 else { 1346 float rry = rx*co + ry*si; 1347 float rrx = -rx*si + ry*co; 1348 dx += rrx; 1349 dy += rry; 1350 } 1351 } 1352 nsamples++; 1353 } 1354 } 1355 di /= nsamples; 1356 dx /= nsamples; 1357 dy /= nsamples; 1358 1359 values[valpos] = di; 1360 if (chan > 1) { 1361 values[valpos + 1] = dx; 1362 } 1363 if (chan > 2) { 1364 values[valpos + 2] = dy; 1365 } 1366 valpos += chan; 1367 } 1368 } 1369 } 1370 1371 void MLDB_Full_Descriptor_Invoker::MLDB_Binary_Comparisons(float* values, unsigned char* desc, 1372 int count, int& dpos) const { 1373 int chan = options_->descriptor_channels; 1374 int* ivalues = (int*) values; 1375 for(int i = 0; i < count * chan; i++) { 1376 ivalues[i] = CV_TOGGLE_FLT(ivalues[i]); 1377 } 1378 1379 for(int pos = 0; pos < chan; pos++) { 1380 for (int i = 0; i < count; i++) { 1381 int ival = ivalues[chan * i + pos]; 1382 for (int j = i + 1; j < count; j++) { 1383 int res = ival > ivalues[chan * j + pos]; 1384 desc[dpos >> 3] |= (res << (dpos & 7)); 1385 dpos++; 1386 } 1387 } 1388 } 1389 } 1390 1391 /* ************************************************************************* */ 1392 /** 1393 * @brief This method computes the descriptor of the provided keypoint given the 1394 * main orientation of the keypoint 1395 * @param kpt Input keypoint 1396 * @param desc Descriptor vector 1397 */ 1398 void MLDB_Full_Descriptor_Invoker::Get_MLDB_Full_Descriptor(const KeyPoint& kpt, unsigned char *desc) const { 1399 1400 const int max_channels = 3; 1401 CV_Assert(options_->descriptor_channels <= max_channels); 1402 float values[16*max_channels]; 1403 const double size_mult[3] = {1, 2.0/3.0, 1.0/2.0}; 1404 1405 float ratio = (float)(1 << kpt.octave); 1406 float scale = (float)fRound(0.5f*kpt.size / ratio); 1407 float xf = kpt.pt.x / ratio; 1408 float yf = kpt.pt.y / ratio; 1409 float co = cos(kpt.angle); 1410 float si = sin(kpt.angle); 1411 int pattern_size = options_->descriptor_pattern_size; 1412 1413 int dpos = 0; 1414 for(int lvl = 0; lvl < 3; lvl++) { 1415 1416 int val_count = (lvl + 2) * (lvl + 2); 1417 int sample_step = static_cast<int>(ceil(pattern_size * size_mult[lvl])); 1418 MLDB_Fill_Values(values, sample_step, kpt.class_id, xf, yf, co, si, scale); 1419 MLDB_Binary_Comparisons(values, desc, val_count, dpos); 1420 } 1421 } 1422 1423 /* ************************************************************************* */ 1424 /** 1425 * @brief This method computes the M-LDB descriptor of the provided keypoint given the 1426 * main orientation of the keypoint. The descriptor is computed based on a subset of 1427 * the bits of the whole descriptor 1428 * @param kpt Input keypoint 1429 * @param desc Descriptor vector 1430 */ 1431 void MLDB_Descriptor_Subset_Invoker::Get_MLDB_Descriptor_Subset(const KeyPoint& kpt, unsigned char *desc) const { 1432 1433 float di = 0.f, dx = 0.f, dy = 0.f; 1434 float rx = 0.f, ry = 0.f; 1435 float sample_x = 0.f, sample_y = 0.f; 1436 int x1 = 0, y1 = 0; 1437 1438 const AKAZEOptions & options = *options_; 1439 const std::vector<TEvolution>& evolution = *evolution_; 1440 1441 // Get the information from the keypoint 1442 float ratio = (float)(1 << kpt.octave); 1443 int scale = fRound(0.5f*kpt.size / ratio); 1444 float angle = kpt.angle; 1445 int level = kpt.class_id; 1446 float yf = kpt.pt.y / ratio; 1447 float xf = kpt.pt.x / ratio; 1448 float co = cos(angle); 1449 float si = sin(angle); 1450 1451 // Allocate memory for the matrix of values 1452 Mat values = Mat_<float>::zeros((4 + 9 + 16)*options.descriptor_channels, 1); 1453 1454 // Sample everything, but only do the comparisons 1455 vector<int> steps(3); 1456 steps.at(0) = options.descriptor_pattern_size; 1457 steps.at(1) = (int)ceil(2.f*options.descriptor_pattern_size / 3.f); 1458 steps.at(2) = options.descriptor_pattern_size / 2; 1459 1460 for (int i = 0; i < descriptorSamples_.rows; i++) { 1461 const int *coords = descriptorSamples_.ptr<int>(i); 1462 int sample_step = steps.at(coords[0]); 1463 di = 0.0f; 1464 dx = 0.0f; 1465 dy = 0.0f; 1466 1467 for (int k = coords[1]; k < coords[1] + sample_step; k++) { 1468 for (int l = coords[2]; l < coords[2] + sample_step; l++) { 1469 1470 // Get the coordinates of the sample point 1471 sample_y = yf + (l*scale*co + k*scale*si); 1472 sample_x = xf + (-l*scale*si + k*scale*co); 1473 1474 y1 = fRound(sample_y); 1475 x1 = fRound(sample_x); 1476 1477 di += *(evolution[level].Lt.ptr<float>(y1)+x1); 1478 1479 if (options.descriptor_channels > 1) { 1480 rx = *(evolution[level].Lx.ptr<float>(y1)+x1); 1481 ry = *(evolution[level].Ly.ptr<float>(y1)+x1); 1482 1483 if (options.descriptor_channels == 2) { 1484 dx += sqrtf(rx*rx + ry*ry); 1485 } 1486 else if (options.descriptor_channels == 3) { 1487 // Get the x and y derivatives on the rotated axis 1488 dx += rx*co + ry*si; 1489 dy += -rx*si + ry*co; 1490 } 1491 } 1492 } 1493 } 1494 1495 *(values.ptr<float>(options.descriptor_channels*i)) = di; 1496 1497 if (options.descriptor_channels == 2) { 1498 *(values.ptr<float>(options.descriptor_channels*i + 1)) = dx; 1499 } 1500 else if (options.descriptor_channels == 3) { 1501 *(values.ptr<float>(options.descriptor_channels*i + 1)) = dx; 1502 *(values.ptr<float>(options.descriptor_channels*i + 2)) = dy; 1503 } 1504 } 1505 1506 // Do the comparisons 1507 const float *vals = values.ptr<float>(0); 1508 const int *comps = descriptorBits_.ptr<int>(0); 1509 1510 for (int i = 0; i<descriptorBits_.rows; i++) { 1511 if (vals[comps[2 * i]] > vals[comps[2 * i + 1]]) { 1512 desc[i / 8] |= (1 << (i % 8)); 1513 } 1514 } 1515 } 1516 1517 /* ************************************************************************* */ 1518 /** 1519 * @brief This method computes the upright (not rotation invariant) M-LDB descriptor 1520 * of the provided keypoint given the main orientation of the keypoint. 1521 * The descriptor is computed based on a subset of the bits of the whole descriptor 1522 * @param kpt Input keypoint 1523 * @param desc Descriptor vector 1524 */ 1525 void Upright_MLDB_Descriptor_Subset_Invoker::Get_Upright_MLDB_Descriptor_Subset(const KeyPoint& kpt, unsigned char *desc) const { 1526 1527 float di = 0.0f, dx = 0.0f, dy = 0.0f; 1528 float rx = 0.0f, ry = 0.0f; 1529 float sample_x = 0.0f, sample_y = 0.0f; 1530 int x1 = 0, y1 = 0; 1531 1532 const AKAZEOptions & options = *options_; 1533 const std::vector<TEvolution>& evolution = *evolution_; 1534 1535 // Get the information from the keypoint 1536 float ratio = (float)(1 << kpt.octave); 1537 int scale = fRound(0.5f*kpt.size / ratio); 1538 int level = kpt.class_id; 1539 float yf = kpt.pt.y / ratio; 1540 float xf = kpt.pt.x / ratio; 1541 1542 // Allocate memory for the matrix of values 1543 Mat values = Mat_<float>::zeros((4 + 9 + 16)*options.descriptor_channels, 1); 1544 1545 vector<int> steps(3); 1546 steps.at(0) = options.descriptor_pattern_size; 1547 steps.at(1) = static_cast<int>(ceil(2.f*options.descriptor_pattern_size / 3.f)); 1548 steps.at(2) = options.descriptor_pattern_size / 2; 1549 1550 for (int i = 0; i < descriptorSamples_.rows; i++) { 1551 const int *coords = descriptorSamples_.ptr<int>(i); 1552 int sample_step = steps.at(coords[0]); 1553 di = 0.0f, dx = 0.0f, dy = 0.0f; 1554 1555 for (int k = coords[1]; k < coords[1] + sample_step; k++) { 1556 for (int l = coords[2]; l < coords[2] + sample_step; l++) { 1557 1558 // Get the coordinates of the sample point 1559 sample_y = yf + l*scale; 1560 sample_x = xf + k*scale; 1561 1562 y1 = fRound(sample_y); 1563 x1 = fRound(sample_x); 1564 di += *(evolution[level].Lt.ptr<float>(y1)+x1); 1565 1566 if (options.descriptor_channels > 1) { 1567 rx = *(evolution[level].Lx.ptr<float>(y1)+x1); 1568 ry = *(evolution[level].Ly.ptr<float>(y1)+x1); 1569 1570 if (options.descriptor_channels == 2) { 1571 dx += sqrtf(rx*rx + ry*ry); 1572 } 1573 else if (options.descriptor_channels == 3) { 1574 dx += rx; 1575 dy += ry; 1576 } 1577 } 1578 } 1579 } 1580 1581 *(values.ptr<float>(options.descriptor_channels*i)) = di; 1582 1583 if (options.descriptor_channels == 2) { 1584 *(values.ptr<float>(options.descriptor_channels*i + 1)) = dx; 1585 } 1586 else if (options.descriptor_channels == 3) { 1587 *(values.ptr<float>(options.descriptor_channels*i + 1)) = dx; 1588 *(values.ptr<float>(options.descriptor_channels*i + 2)) = dy; 1589 } 1590 } 1591 1592 // Do the comparisons 1593 const float *vals = values.ptr<float>(0); 1594 const int *comps = descriptorBits_.ptr<int>(0); 1595 1596 for (int i = 0; i<descriptorBits_.rows; i++) { 1597 if (vals[comps[2 * i]] > vals[comps[2 * i + 1]]) { 1598 desc[i / 8] |= (1 << (i % 8)); 1599 } 1600 } 1601 } 1602 1603 /* ************************************************************************* */ 1604 /** 1605 * @brief This function computes a (quasi-random) list of bits to be taken 1606 * from the full descriptor. To speed the extraction, the function creates 1607 * a list of the samples that are involved in generating at least a bit (sampleList) 1608 * and a list of the comparisons between those samples (comparisons) 1609 * @param sampleList 1610 * @param comparisons The matrix with the binary comparisons 1611 * @param nbits The number of bits of the descriptor 1612 * @param pattern_size The pattern size for the binary descriptor 1613 * @param nchannels Number of channels to consider in the descriptor (1-3) 1614 * @note The function keeps the 18 bits (3-channels by 6 comparisons) of the 1615 * coarser grid, since it provides the most robust estimations 1616 */ 1617 void generateDescriptorSubsample(Mat& sampleList, Mat& comparisons, int nbits, 1618 int pattern_size, int nchannels) { 1619 1620 int ssz = 0; 1621 for (int i = 0; i < 3; i++) { 1622 int gz = (i + 2)*(i + 2); 1623 ssz += gz*(gz - 1) / 2; 1624 } 1625 ssz *= nchannels; 1626 1627 CV_Assert(nbits <= ssz); // Descriptor size can't be bigger than full descriptor 1628 1629 // Since the full descriptor is usually under 10k elements, we pick 1630 // the selection from the full matrix. We take as many samples per 1631 // pick as the number of channels. For every pick, we 1632 // take the two samples involved and put them in the sampling list 1633 1634 Mat_<int> fullM(ssz / nchannels, 5); 1635 for (int i = 0, c = 0; i < 3; i++) { 1636 int gdiv = i + 2; //grid divisions, per row 1637 int gsz = gdiv*gdiv; 1638 int psz = (int)ceil(2.f*pattern_size / (float)gdiv); 1639 1640 for (int j = 0; j < gsz; j++) { 1641 for (int k = j + 1; k < gsz; k++, c++) { 1642 fullM(c, 0) = i; 1643 fullM(c, 1) = psz*(j % gdiv) - pattern_size; 1644 fullM(c, 2) = psz*(j / gdiv) - pattern_size; 1645 fullM(c, 3) = psz*(k % gdiv) - pattern_size; 1646 fullM(c, 4) = psz*(k / gdiv) - pattern_size; 1647 } 1648 } 1649 } 1650 1651 srand(1024); 1652 Mat_<int> comps = Mat_<int>(nchannels * (int)ceil(nbits / (float)nchannels), 2); 1653 comps = 1000; 1654 1655 // Select some samples. A sample includes all channels 1656 int count = 0; 1657 int npicks = (int)ceil(nbits / (float)nchannels); 1658 Mat_<int> samples(29, 3); 1659 Mat_<int> fullcopy = fullM.clone(); 1660 samples = -1; 1661 1662 for (int i = 0; i < npicks; i++) { 1663 int k = rand() % (fullM.rows - i); 1664 if (i < 6) { 1665 // Force use of the coarser grid values and comparisons 1666 k = i; 1667 } 1668 1669 bool n = true; 1670 1671 for (int j = 0; j < count; j++) { 1672 if (samples(j, 0) == fullcopy(k, 0) && samples(j, 1) == fullcopy(k, 1) && samples(j, 2) == fullcopy(k, 2)) { 1673 n = false; 1674 comps(i*nchannels, 0) = nchannels*j; 1675 comps(i*nchannels + 1, 0) = nchannels*j + 1; 1676 comps(i*nchannels + 2, 0) = nchannels*j + 2; 1677 break; 1678 } 1679 } 1680 1681 if (n) { 1682 samples(count, 0) = fullcopy(k, 0); 1683 samples(count, 1) = fullcopy(k, 1); 1684 samples(count, 2) = fullcopy(k, 2); 1685 comps(i*nchannels, 0) = nchannels*count; 1686 comps(i*nchannels + 1, 0) = nchannels*count + 1; 1687 comps(i*nchannels + 2, 0) = nchannels*count + 2; 1688 count++; 1689 } 1690 1691 n = true; 1692 for (int j = 0; j < count; j++) { 1693 if (samples(j, 0) == fullcopy(k, 0) && samples(j, 1) == fullcopy(k, 3) && samples(j, 2) == fullcopy(k, 4)) { 1694 n = false; 1695 comps(i*nchannels, 1) = nchannels*j; 1696 comps(i*nchannels + 1, 1) = nchannels*j + 1; 1697 comps(i*nchannels + 2, 1) = nchannels*j + 2; 1698 break; 1699 } 1700 } 1701 1702 if (n) { 1703 samples(count, 0) = fullcopy(k, 0); 1704 samples(count, 1) = fullcopy(k, 3); 1705 samples(count, 2) = fullcopy(k, 4); 1706 comps(i*nchannels, 1) = nchannels*count; 1707 comps(i*nchannels + 1, 1) = nchannels*count + 1; 1708 comps(i*nchannels + 2, 1) = nchannels*count + 2; 1709 count++; 1710 } 1711 1712 Mat tmp = fullcopy.row(k); 1713 fullcopy.row(fullcopy.rows - i - 1).copyTo(tmp); 1714 } 1715 1716 sampleList = samples.rowRange(0, count).clone(); 1717 comparisons = comps.rowRange(0, nbits).clone(); 1718 } 1719 1720 } 1721
