From 7a594354903f7a44239854df2763d3c8f1505b56 Mon Sep 17 00:00:00 2001 From: Ievgen Khvedchenia Date: Fri, 4 Apr 2014 14:25:38 +0300 Subject: [PATCH] KAZE and AKAZE integration initial commit --- modules/features2d/src/akaze/AKAZE.cpp | 2068 ++++++++++++ modules/features2d/src/akaze/AKAZE.h | 175 + modules/features2d/src/akaze/config.h | 155 + modules/features2d/src/akaze/fed.h | 26 + .../src/akaze/nldiffusion_functions.cpp | 431 +++ .../src/akaze/nldiffusion_functions.h | 41 + modules/features2d/src/akaze/utils.cpp | 196 ++ modules/features2d/src/akaze/utils.h | 54 + modules/features2d/src/kaze.cpp | 0 modules/features2d/src/kaze/KAZE.cpp | 2801 +++++++++++++++++ modules/features2d/src/kaze/KAZE.h | 294 ++ modules/features2d/src/kaze/config.h | 129 + modules/features2d/src/kaze/fed.cpp | 192 ++ modules/features2d/src/kaze/fed.h | 30 + .../src/kaze/nldiffusion_functions.cpp | 386 +++ .../src/kaze/nldiffusion_functions.h | 51 + modules/features2d/src/kaze/utils.cpp | 92 + modules/features2d/src/kaze/utils.h | 41 + 18 files changed, 7162 insertions(+) create mode 100644 modules/features2d/src/akaze/AKAZE.cpp create mode 100644 modules/features2d/src/akaze/AKAZE.h create mode 100644 modules/features2d/src/akaze/config.h create mode 100644 modules/features2d/src/akaze/fed.h create mode 100644 modules/features2d/src/akaze/nldiffusion_functions.cpp create mode 100644 modules/features2d/src/akaze/nldiffusion_functions.h create mode 100644 modules/features2d/src/akaze/utils.cpp create mode 100644 modules/features2d/src/akaze/utils.h create mode 100644 modules/features2d/src/kaze.cpp create mode 100644 modules/features2d/src/kaze/KAZE.cpp create mode 100755 modules/features2d/src/kaze/KAZE.h create mode 100644 modules/features2d/src/kaze/config.h create mode 100644 modules/features2d/src/kaze/fed.cpp create mode 100644 modules/features2d/src/kaze/fed.h create mode 100644 modules/features2d/src/kaze/nldiffusion_functions.cpp create mode 100755 modules/features2d/src/kaze/nldiffusion_functions.h create mode 100644 modules/features2d/src/kaze/utils.cpp create mode 100644 modules/features2d/src/kaze/utils.h diff --git a/modules/features2d/src/akaze/AKAZE.cpp b/modules/features2d/src/akaze/AKAZE.cpp new file mode 100644 index 0000000000..5a110ac175 --- /dev/null +++ b/modules/features2d/src/akaze/AKAZE.cpp @@ -0,0 +1,2068 @@ +//============================================================================= +// +// AKAZE.cpp +// Authors: Pablo F. Alcantarilla (1), Jesus Nuevo (2) +// Institutions: Georgia Institute of Technology (1) +// TrueVision Solutions (2) +// Date: 15/09/2013 +// Email: pablofdezalc@gmail.com +// +// AKAZE Features Copyright 2013, Pablo F. Alcantarilla, Jesus Nuevo +// All Rights Reserved +// See LICENSE for the license information +//============================================================================= + +/** + * @file AKAZE.cpp + * @brief Main class for detecting and describing binary features in an + * accelerated nonlinear scale space + * @date Sep 15, 2013 + * @author Pablo F. Alcantarilla, Jesus Nuevo + */ + +#include "AKAZE.h" + +using namespace std; +using namespace cv; + +//******************************************************************************* +//******************************************************************************* + +/** + * @brief AKAZE constructor with input options + * @param options AKAZE configuration options + * @note This constructor allocates memory for the nonlinear scale space +*/ +AKAZE::AKAZE(const AKAZEOptions& options) { + + soffset_ = options.soffset; + factor_size_ = DEFAULT_FACTOR_SIZE; + sderivatives_ = DEFAULT_SIGMA_SMOOTHING_DERIVATIVES; + omax_ = options.omax; + nsublevels_ = options.nsublevels; + dthreshold_ = options.dthreshold; + descriptor_ = options.descriptor; + diffusivity_ = options.diffusivity; + save_scale_space_ = options.save_scale_space; + verbosity_ = options.verbosity; + img_width_ = options.img_width; + img_height_ = options.img_height; + noctaves_ = omax_; + ncycles_ = 0; + reordering_ = true; + descriptor_size_ = options.descriptor_size; + descriptor_channels_ = options.descriptor_channels; + descriptor_pattern_size_ = options.descriptor_pattern_size; + tkcontrast_ = 0.0; + tscale_ = 0.0; + tderivatives_ = 0.0; + tdetector_ = 0.0; + textrema_ = 0.0; + tsubpixel_ = 0.0; + tdescriptor_ = 0.0; + + if (descriptor_size_ > 0 && descriptor_ >= MLDB_UPRIGHT) { + generateDescriptorSubsample(descriptorSamples_,descriptorBits_,descriptor_size_, + descriptor_pattern_size_,descriptor_channels_); + } + + Allocate_Memory_Evolution(); +} + +//******************************************************************************* +//******************************************************************************* + +/** + * @brief AKAZE destructor +*/ +AKAZE::~AKAZE(void) { + + evolution_.clear(); +} + +//******************************************************************************* +//******************************************************************************* + +/** + * @brief This method allocates the memory for the nonlinear diffusion evolution +*/ +void AKAZE::Allocate_Memory_Evolution(void) { + + float rfactor = 0.0; + int level_height = 0, level_width = 0; + + // Allocate the dimension of the matrices for the evolution + for (int i = 0; i <= omax_-1 && i <= DEFAULT_OCTAVE_MAX; i++) { + rfactor = 1.0/pow(2.f,i); + level_height = (int)(img_height_*rfactor); + level_width = (int)(img_width_*rfactor); + + // Smallest possible octave + if (level_width < 80 || level_height < 40) { + noctaves_ = i; + i = omax_; + break; + } + + for (int j = 0; j < nsublevels_; j++) { + tevolution aux; + aux.Lx = cv::Mat::zeros(level_height,level_width,CV_32F); + aux.Ly = cv::Mat::zeros(level_height,level_width,CV_32F); + aux.Lxx = cv::Mat::zeros(level_height,level_width,CV_32F); + aux.Lxy = cv::Mat::zeros(level_height,level_width,CV_32F); + aux.Lyy = cv::Mat::zeros(level_height,level_width,CV_32F); + aux.Lt = cv::Mat::zeros(level_height,level_width,CV_32F); + aux.Ldet = cv::Mat::zeros(level_height,level_width,CV_32F); + aux.Lflow = cv::Mat::zeros(level_height,level_width,CV_32F); + aux.Lstep = cv::Mat::zeros(level_height,level_width,CV_32F); + aux.esigma = soffset_*pow(2.f,(float)(j)/(float)(nsublevels_) + i); + aux.sigma_size = fRound(aux.esigma); + aux.etime = 0.5*(aux.esigma*aux.esigma); + aux.octave = i; + aux.sublevel = j; + evolution_.push_back(aux); + } + } + + // Allocate memory for the number of cycles and time steps + for (size_t i = 1; i < evolution_.size(); i++) { + int naux = 0; + std::vector tau; + float ttime = 0.0; + ttime = evolution_[i].etime-evolution_[i-1].etime; + naux = fed_tau_by_process_time(ttime,1,0.25,reordering_,tau); + nsteps_.push_back(naux); + tsteps_.push_back(tau); + ncycles_++; + } +} + +//******************************************************************************* +//******************************************************************************* + +/** + * @brief This method creates the nonlinear scale space for a given image + * @param img Input image for which the nonlinear scale space needs to be created + * @return 0 if the nonlinear scale space was created successfully, -1 otherwise +*/ +int AKAZE::Create_Nonlinear_Scale_Space(const cv::Mat &img) { + + double t1 = 0.0, t2 = 0.0; + + if (evolution_.size() == 0) { + cout << "Error generating the nonlinear scale space!!" << endl; + cout << "Firstly you need to call AKAZE::Allocate_Memory_Evolution()" << endl; + return -1; + } + + t1 = getTickCount(); + + // Copy the original image to the first level of the evolution + img.copyTo(evolution_[0].Lt); + gaussian_2D_convolution(evolution_[0].Lt,evolution_[0].Lt,0,0,soffset_); + evolution_[0].Lt.copyTo(evolution_[0].Lsmooth); + + // Firstly compute the kcontrast factor + kcontrast_ = compute_k_percentile(img,KCONTRAST_PERCENTILE,1.0,KCONTRAST_NBINS,0,0); + + t2 = getTickCount(); + tkcontrast_ = 1000.0*(t2-t1) / getTickFrequency(); + + // Now generate the rest of evolution levels + for (size_t i = 1; i < evolution_.size(); i++) { + + if (evolution_[i].octave > evolution_[i-1].octave) { + halfsample_image(evolution_[i-1].Lt,evolution_[i].Lt); + kcontrast_ = kcontrast_*0.75; + } + else { + evolution_[i-1].Lt.copyTo(evolution_[i].Lt); + } + + gaussian_2D_convolution(evolution_[i].Lt,evolution_[i].Lsmooth,0,0,1.0); + + // Compute the Gaussian derivatives Lx and Ly + image_derivatives_scharr(evolution_[i].Lsmooth,evolution_[i].Lx,1,0); + image_derivatives_scharr(evolution_[i].Lsmooth,evolution_[i].Ly,0,1); + + // Compute the conductivity equation + switch (diffusivity_) { + case 0: + pm_g1(evolution_[i].Lx,evolution_[i].Ly,evolution_[i].Lflow,kcontrast_); + break; + case 1: + pm_g2(evolution_[i].Lx,evolution_[i].Ly,evolution_[i].Lflow,kcontrast_); + break; + case 2: + weickert_diffusivity(evolution_[i].Lx,evolution_[i].Ly,evolution_[i].Lflow,kcontrast_); + break; + case 3: + charbonnier_diffusivity(evolution_[i].Lx,evolution_[i].Ly,evolution_[i].Lflow,kcontrast_); + break; + default: + std::cerr << "Diffusivity: " << diffusivity_ << " is not supported" << std::endl; + } + + // Perform FED n inner steps + for (int j = 0; j < nsteps_[i-1]; j++) { + nld_step_scalar(evolution_[i].Lt,evolution_[i].Lflow,evolution_[i].Lstep,tsteps_[i-1][j]); + } + } + + t2 = getTickCount(); + tscale_ = 1000.0*(t2-t1) / getTickFrequency(); + + return 0; +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method selects interesting keypoints through the nonlinear scale space + * @param kpts Vector of detected keypoints +*/ +void AKAZE::Feature_Detection(std::vector& kpts) { + + double t1 = 0.0, t2 = 0.0; + + t1 = getTickCount(); + + Compute_Determinant_Hessian_Response(); + Find_Scale_Space_Extrema(kpts); + Do_Subpixel_Refinement(kpts); + + t2 = getTickCount(); + tdetector_ = 1000.0*(t2-t1) / getTickFrequency(); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the multiscale derivatives for the nonlinear scale space +*/ +void AKAZE::Compute_Multiscale_Derivatives(void) { + + double t1 = 0.0, t2 = 0.0; + + t1 = getTickCount(); + +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (size_t i = 0; i < evolution_.size(); i++) { + float ratio = pow(2.f,evolution_[i].octave); + int sigma_size_ = fRound(evolution_[i].esigma*factor_size_/ratio); + + compute_scharr_derivatives(evolution_[i].Lsmooth,evolution_[i].Lx,1,0,sigma_size_); + compute_scharr_derivatives(evolution_[i].Lsmooth,evolution_[i].Ly,0,1,sigma_size_); + compute_scharr_derivatives(evolution_[i].Lx,evolution_[i].Lxx,1,0,sigma_size_); + compute_scharr_derivatives(evolution_[i].Ly,evolution_[i].Lyy,0,1,sigma_size_); + compute_scharr_derivatives(evolution_[i].Lx,evolution_[i].Lxy,0,1,sigma_size_); + + evolution_[i].Lx = evolution_[i].Lx*((sigma_size_)); + evolution_[i].Ly = evolution_[i].Ly*((sigma_size_)); + evolution_[i].Lxx = evolution_[i].Lxx*((sigma_size_)*(sigma_size_)); + evolution_[i].Lxy = evolution_[i].Lxy*((sigma_size_)*(sigma_size_)); + evolution_[i].Lyy = evolution_[i].Lyy*((sigma_size_)*(sigma_size_)); + } + + t2 = getTickCount(); + tderivatives_ = 1000.0*(t2-t1) / getTickFrequency(); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the feature detector response for the nonlinear scale space + * @note We use the Hessian determinant as the feature detector response +*/ +void AKAZE::Compute_Determinant_Hessian_Response(void) { + + // Firstly compute the multiscale derivatives + Compute_Multiscale_Derivatives(); + + for (size_t i = 0; i < evolution_.size(); i++) { + if (verbosity_ == true) { + cout << "Computing detector response. Determinant of Hessian. Evolution time: " << evolution_[i].etime << endl; + } + + for (int ix = 0; ix < evolution_[i].Ldet.rows; ix++) { + for (int jx = 0; jx < evolution_[i].Ldet.cols; jx++) { + float lxx = *(evolution_[i].Lxx.ptr(ix)+jx); + float lxy = *(evolution_[i].Lxy.ptr(ix)+jx); + float lyy = *(evolution_[i].Lyy.ptr(ix)+jx); + *(evolution_[i].Ldet.ptr(ix)+jx) = (lxx*lyy-lxy*lxy); + } + } + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method finds extrema in the nonlinear scale space + * @param kpts Vector of detected keypoints +*/ +void AKAZE::Find_Scale_Space_Extrema(std::vector& kpts) { + + double t1 = 0.0, t2 = 0.0; + float value = 0.0; + float dist = 0.0, ratio = 0.0, smax = 0.0; + int npoints = 0, id_repeated = 0; + int sigma_size_ = 0, left_x = 0, right_x = 0, up_y = 0, down_y = 0; + bool is_extremum = false, is_repeated = false, is_out = false; + cv::KeyPoint point; + + // Set maximum size + if (descriptor_ == SURF_UPRIGHT || descriptor_ == SURF || + descriptor_ == MLDB_UPRIGHT || descriptor_ == MLDB) { + smax = 10.0*sqrtf(2.0); + } + else if (descriptor_ == MSURF_UPRIGHT || descriptor_ == MSURF) { + smax = 12.0*sqrtf(2.0); + } + + t1 = getTickCount(); + + for (size_t i = 0; i < evolution_.size(); i++) { + for (int ix = 1; ix < evolution_[i].Ldet.rows-1; ix++) { + for (int jx = 1; jx < evolution_[i].Ldet.cols-1; jx++) { + is_extremum = false; + is_repeated = false; + is_out = false; + value = *(evolution_[i].Ldet.ptr(ix)+jx); + + // Filter the points with the detector threshold + if (value > dthreshold_ && value >= DEFAULT_MIN_DETECTOR_THRESHOLD && + value > *(evolution_[i].Ldet.ptr(ix)+jx-1) && + value > *(evolution_[i].Ldet.ptr(ix)+jx+1) && + value > *(evolution_[i].Ldet.ptr(ix-1)+jx-1) && + value > *(evolution_[i].Ldet.ptr(ix-1)+jx) && + value > *(evolution_[i].Ldet.ptr(ix-1)+jx+1) && + value > *(evolution_[i].Ldet.ptr(ix+1)+jx-1) && + value > *(evolution_[i].Ldet.ptr(ix+1)+jx) && + value > *(evolution_[i].Ldet.ptr(ix+1)+jx+1)) { + is_extremum = true; + + point.response = fabs(value); + point.size = evolution_[i].esigma*factor_size_; + point.octave = evolution_[i].octave; + point.class_id = i; + ratio = pow(2.f,point.octave); + sigma_size_ = fRound(point.size/ratio); + point.pt.x = jx; + point.pt.y = ix; + + for (size_t ik = 0; ik < kpts.size(); ik++) { + if (point.class_id == kpts[ik].class_id-1 || + point.class_id == kpts[ik].class_id || + point.class_id == kpts[ik].class_id+1) { + dist = sqrt(pow(point.pt.x*ratio-kpts[ik].pt.x,2)+pow(point.pt.y*ratio-kpts[ik].pt.y,2)); + if (dist <= point.size) { + if (point.response > kpts[ik].response) { + id_repeated = ik; + is_repeated = true; + } + else { + is_extremum = false; + } + break; + } + } + } + + // Check out of bounds + if (is_extremum == true) { + // Check that the point is under the image limits for the descriptor computation + left_x = fRound(point.pt.x-smax*sigma_size_)-1; + right_x = fRound(point.pt.x+smax*sigma_size_) +1; + up_y = fRound(point.pt.y-smax*sigma_size_)-1; + down_y = fRound(point.pt.y+smax*sigma_size_)+1; + + if (left_x < 0 || right_x >= evolution_[i].Ldet.cols || + up_y < 0 || down_y >= evolution_[i].Ldet.rows) { + is_out = true; + } + + if (is_out == false) { + if (is_repeated == false) { + point.pt.x *= ratio; + point.pt.y *= ratio; + kpts.push_back(point); + npoints++; + } + else { + point.pt.x *= ratio; + point.pt.y *= ratio; + kpts[id_repeated] = point; + } + } // if is_out + } //if is_extremum + } + } // for jx + } // for ix + } // for i + + t2 = getTickCount(); + textrema_ = 1000.0*(t2-t1) / getTickFrequency(); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method performs subpixel refinement of the detected keypoints + * @param kpts Vector of detected keypoints +*/ +void AKAZE::Do_Subpixel_Refinement(std::vector& kpts) { + + double t1 = 0.0, t2 = 0.0; + float Dx = 0.0, Dy = 0.0, ratio = 0.0; + float Dxx = 0.0, Dyy = 0.0, Dxy = 0.0; + int x = 0, y = 0; + Mat A = Mat::zeros(2,2,CV_32F); + Mat b = Mat::zeros(2,1,CV_32F); + Mat dst = Mat::zeros(2,1,CV_32F); + + t1 = getTickCount(); + + for (size_t i = 0; i < kpts.size(); i++) { + ratio = pow(2.f,kpts[i].octave); + x = fRound(kpts[i].pt.x/ratio); + y = fRound(kpts[i].pt.y/ratio); + + // Compute the gradient + Dx = (0.5)*(*(evolution_[kpts[i].class_id].Ldet.ptr(y)+x+1) + -*(evolution_[kpts[i].class_id].Ldet.ptr(y)+x-1)); + Dy = (0.5)*(*(evolution_[kpts[i].class_id].Ldet.ptr(y+1)+x) + -*(evolution_[kpts[i].class_id].Ldet.ptr(y-1)+x)); + + // Compute the Hessian + Dxx = (*(evolution_[kpts[i].class_id].Ldet.ptr(y)+x+1) + + *(evolution_[kpts[i].class_id].Ldet.ptr(y)+x-1) + -2.0*(*(evolution_[kpts[i].class_id].Ldet.ptr(y)+x))); + + Dyy = (*(evolution_[kpts[i].class_id].Ldet.ptr(y+1)+x) + + *(evolution_[kpts[i].class_id].Ldet.ptr(y-1)+x) + -2.0*(*(evolution_[kpts[i].class_id].Ldet.ptr(y)+x))); + + Dxy = (0.25)*(*(evolution_[kpts[i].class_id].Ldet.ptr(y+1)+x+1) + +(*(evolution_[kpts[i].class_id].Ldet.ptr(y-1)+x-1))) + -(0.25)*(*(evolution_[kpts[i].class_id].Ldet.ptr(y-1)+x+1) + +(*(evolution_[kpts[i].class_id].Ldet.ptr(y+1)+x-1))); + + // Solve the linear system + *(A.ptr(0)) = Dxx; + *(A.ptr(1)+1) = Dyy; + *(A.ptr(0)+1) = *(A.ptr(1)) = Dxy; + *(b.ptr(0)) = -Dx; + *(b.ptr(1)) = -Dy; + + solve(A,b,dst,DECOMP_LU); + + if (fabs(*(dst.ptr(0))) <= 1.0 && fabs(*(dst.ptr(1))) <= 1.0) { + kpts[i].pt.x = x + (*(dst.ptr(0))); + kpts[i].pt.y = y + (*(dst.ptr(1))); + kpts[i].pt.x *= powf(2.f,evolution_[kpts[i].class_id].octave); + kpts[i].pt.y *= powf(2.f,evolution_[kpts[i].class_id].octave); + kpts[i].angle = 0.0; + + // In OpenCV the size of a keypoint its the diameter + kpts[i].size *= 2.0; + } + // Delete the point since its not stable + else { + kpts.erase(kpts.begin()+i); + i--; + } + } + + t2 = getTickCount(); + tsubpixel_ = 1000.0*(t2-t1) / getTickFrequency(); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method performs feature suppression based on 2D distance + * @param kpts Vector of keypoints + * @param mdist Maximum distance in pixels +*/ +void AKAZE::Feature_Suppression_Distance(std::vector& kpts, float mdist) { + + vector aux; + vector to_delete; + float dist = 0.0, x1 = 0.0, y1 = 0.0, x2 = 0.0, y2 = 0.0; + bool found = false; + + for (size_t i = 0; i < kpts.size(); i++) { + x1 = kpts[i].pt.x; + y1 = kpts[i].pt.y; + for (size_t j = i+1; j < kpts.size(); j++) { + x2 = kpts[j].pt.x; + y2 = kpts[j].pt.y; + dist = sqrt(pow(x1-x2,2)+pow(y1-y2,2)); + if (dist < mdist) { + if (fabs(kpts[i].response) >= fabs(kpts[j].response)) { + to_delete.push_back(j); + } + else { + to_delete.push_back(i); + break; + } + } + } + } + + for (size_t i = 0; i < kpts.size(); i++) { + found = false; + for (size_t j = 0; j < to_delete.size(); j++) { + if (i == (size_t)(to_delete[j])) { + found = true; + break; + } + } + if (found == false) { + aux.push_back(kpts[i]); + } + } + + kpts.clear(); + kpts = aux; + aux.clear(); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the set of descriptors through the nonlinear scale space + * @param kpts Vector of detected keypoints + * @param desc Matrix to store the descriptors +*/ +void AKAZE::Compute_Descriptors(std::vector& kpts, cv::Mat& desc) { + + double t1 = 0.0, t2 = 0.0; + + t1 = getTickCount(); + + // Allocate memory for the matrix with the descriptors + if (descriptor_ < MLDB_UPRIGHT) { + desc = cv::Mat::zeros(kpts.size(),64,CV_32FC1); + } + else { + // We use the full length binary descriptor -> 486 bits + if (descriptor_size_ == 0) { + int t = (6+36+120)*descriptor_channels_; + desc = cv::Mat::zeros(kpts.size(),ceil(t/8.),CV_8UC1); + } + else { + // We use the random bit selection length binary descriptor + desc = cv::Mat::zeros(kpts.size(),ceil(descriptor_size_/8.),CV_8UC1); + } + } + + switch (descriptor_) + { + case SURF_UPRIGHT : // Upright descriptors, not invariant to rotation + { +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (size_t i = 0; i < kpts.size(); i++) { + Get_SURF_Descriptor_Upright_64(kpts[i],desc.ptr(i)); + } + } + break; + case SURF : + { +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (size_t i = 0; i < kpts.size(); i++) { + Compute_Main_Orientation_SURF(kpts[i]); + Get_SURF_Descriptor_64(kpts[i],desc.ptr(i)); + } + } + break; + case MSURF_UPRIGHT : // Upright descriptors, not invariant to rotation + { +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (size_t i = 0; i < kpts.size(); i++) { + Get_MSURF_Upright_Descriptor_64(kpts[i],desc.ptr(i)); + } + } + break; + case MSURF : + { +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (size_t i = 0; i < kpts.size(); i++) { + Compute_Main_Orientation_SURF(kpts[i]); + Get_MSURF_Descriptor_64(kpts[i],desc.ptr(i)); + } + } + break; + case MLDB_UPRIGHT : // Upright descriptors, not invariant to rotation + { +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (size_t i = 0; i < kpts.size(); i++) { + if (descriptor_size_ == 0) + Get_Upright_MLDB_Full_Descriptor(kpts[i],desc.ptr(i)); + else + Get_Upright_MLDB_Descriptor_Subset(kpts[i],desc.ptr(i)); + } + } + break; + case MLDB : + { +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (size_t i = 0; i < kpts.size(); i++) { + Compute_Main_Orientation_SURF(kpts[i]); + if (descriptor_size_ == 0) + Get_MLDB_Full_Descriptor(kpts[i],desc.ptr(i)); + else + Get_MLDB_Descriptor_Subset(kpts[i],desc.ptr(i)); + } + } + break; + } + + t2 = getTickCount(); + tdescriptor_ = 1000.0*(t2-t1) / getTickFrequency(); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the main orientation for a given keypoint + * @param kpt Input keypoint + * @note The orientation is computed using a similar approach as described in the + * original SURF method. See Bay et al., Speeded Up Robust Features, ECCV 2006 +*/ +void AKAZE::Compute_Main_Orientation_SURF(cv::KeyPoint& kpt) { + + int ix = 0, iy = 0, idx = 0, s = 0, level = 0; + float xf = 0.0, yf = 0.0, gweight = 0.0, ratio = 0.0; + std::vector resX(109), resY(109), Ang(109); + const int id[] = {6,5,4,3,2,1,0,1,2,3,4,5,6}; + + // Variables for computing the dominant direction + float sumX = 0.0, sumY = 0.0, max = 0.0, ang1 = 0.0, ang2 = 0.0; + + // Get the information from the keypoint + level = kpt.class_id; + ratio = (float)(1<(iy)+ix)); + resY[idx] = gweight*(*(evolution_[level].Ly.ptr(iy)+ix)); + + Ang[idx] = get_angle(resX[idx],resY[idx]); + ++idx; + } + } + } + + // Loop slides pi/3 window around feature point + for (ang1 = 0; ang1 < 2.0*CV_PI; ang1+=0.15f) { + ang2 =(ang1+CV_PI/3.0f > 2.0*CV_PI ? ang1-5.0f*CV_PI/3.0f : ang1+CV_PI/3.0f); + sumX = sumY = 0.f; + + for (size_t k = 0; k < Ang.size(); ++k) { + // Get angle from the x-axis of the sample point + const float & ang = Ang[k]; + + // Determine whether the point is within the window + if (ang1 < ang2 && ang1 < ang && ang < ang2) { + sumX+=resX[k]; + sumY+=resY[k]; + } + else if (ang2 < ang1 && + ((ang > 0 && ang < ang2) || (ang > ang1 && ang < 2.0*CV_PI) )) { + sumX+=resX[k]; + sumY+=resY[k]; + } + } + + // if the vector produced from this window is longer than all + // previous vectors then this forms the new dominant direction + if (sumX*sumX + sumY*sumY > max) { + // store largest orientation + max = sumX*sumX + sumY*sumY; + kpt.angle = get_angle(sumX, sumY); + } + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the upright descriptor of the provided keypoint + * @param kpt Input keypoint + * @note Rectangular grid of 20 s x 20 s. Descriptor Length 64. No additional + * Gaussian weighting is performed. The descriptor is inspired from Bay et al., + * Speeded Up Robust Features, ECCV, 2006 +*/ +void AKAZE::Get_SURF_Descriptor_Upright_64(const cv::KeyPoint& kpt, float *desc) { + + float dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0; + float rx = 0.0, ry = 0.0, len = 0.0, xf = 0.0, yf = 0.0; + float sample_x = 0.0, sample_y = 0.0; + float fx = 0.0, fy = 0.0, ratio = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0, res4 = 0.0; + int x1 = 0, y1 = 0, x2 = 0, y2 = 0, sample_step = 0, pattern_size = 0, dcount = 0; + int scale = 0, dsize = 0, level = 0; + + // Set the descriptor size and the sample and pattern sizes + dsize = 64; + sample_step = 5; + pattern_size = 10; + + // Get the information from the keypoint + ratio = (float)(1<(y1)+x1); + res2 = *(evolution_[level].Lx.ptr(y1)+x2); + res3 = *(evolution_[level].Lx.ptr(y2)+x1); + res4 = *(evolution_[level].Lx.ptr(y2)+x2); + rx = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Ly.ptr(y1)+x1); + res2 = *(evolution_[level].Ly.ptr(y1)+x2); + res3 = *(evolution_[level].Ly.ptr(y2)+x1); + res4 = *(evolution_[level].Ly.ptr(y2)+x2); + ry = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + // Sum the derivatives to the cumulative descriptor + dx += rx; + dy += ry; + mdx += fabs(rx); + mdy += fabs(ry); + } + } + + // Add the values to the descriptor vector + desc[dcount++] = dx; + desc[dcount++] = dy; + desc[dcount++] = mdx; + desc[dcount++] = mdy; + + // Store the current length^2 of the vector + len += dx*dx + dy*dy + mdx*mdx + mdy*mdy; + } + } + + // convert to unit vector + len = sqrt(len); + + for (int i = 0; i < dsize; i++) { + desc[i] /= len; + } +} + +//************************************************************************************* +//************************************************************************************* +/** + * @brief This method computes the descriptor of the provided keypoint given the + * main orientation + * @param kpt Input keypoint + * @param desc Descriptor vector + * @note Rectangular grid of 20 s x 20 s. Descriptor Length 64. No additional + * Gaussian weighting is performed. The descriptor is inspired from Bay et al., + * Speeded Up Robust Features, ECCV, 2006 +*/ +void AKAZE::Get_SURF_Descriptor_64(const cv::KeyPoint& kpt, float *desc) { + + float dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0; + float rx = 0.0, ry = 0.0, rrx = 0.0, rry = 0.0, len = 0.0, xf = 0.0, yf = 0.0; + float sample_x = 0.0, sample_y = 0.0, co = 0.0, si = 0.0, angle = 0.0; + float fx = 0.0, fy = 0.0, ratio = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0, res4 = 0.0; + int x1 = 0, y1 = 0, x2 = 0, y2 = 0, sample_step = 0, pattern_size = 0, dcount = 0; + int scale = 0, dsize = 0, level = 0; + + // Set the descriptor size and the sample and pattern sizes + dsize = 64; + sample_step = 5; + pattern_size = 10; + + // Get the information from the keypoint + ratio = (float)(1<(y1)+x1); + res2 = *(evolution_[level].Lx.ptr(y1)+x2); + res3 = *(evolution_[level].Lx.ptr(y2)+x1); + res4 = *(evolution_[level].Lx.ptr(y2)+x2); + rx = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Ly.ptr(y1)+x1); + res2 = *(evolution_[level].Ly.ptr(y1)+x2); + res3 = *(evolution_[level].Ly.ptr(y2)+x1); + res4 = *(evolution_[level].Ly.ptr(y2)+x2); + ry = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + // Get the x and y derivatives on the rotated axis + rry = rx*co + ry*si; + rrx = -rx*si + ry*co; + + // Sum the derivatives to the cumulative descriptor + dx += rrx; + dy += rry; + mdx += fabs(rrx); + mdy += fabs(rry); + } + } + + // Add the values to the descriptor vector + desc[dcount++] = dx; + desc[dcount++] = dy; + desc[dcount++] = mdx; + desc[dcount++] = mdy; + + // Store the current length^2 of the vector + len += dx*dx + dy*dy + mdx*mdx + mdy*mdy; + } + } + + // convert to unit vector + len = sqrt(len); + + for (int i = 0; i < dsize; i++) { + desc[i] /= len; + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the upright descriptor (not rotation invariant) of + * the provided keypoint + * @param kpt Input keypoint + * @param desc Descriptor vector + * @note Rectangular grid of 24 s x 24 s. Descriptor Length 64. The descriptor is inspired + * from Agrawal et al., CenSurE: Center Surround Extremas for Realtime Feature Detection and Matching, + * ECCV 2008 +*/ +void AKAZE::Get_MSURF_Upright_Descriptor_64(const cv::KeyPoint& kpt, float *desc) { + + float dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0, gauss_s1 = 0.0, gauss_s2 = 0.0; + float rx = 0.0, ry = 0.0, len = 0.0, xf = 0.0, yf = 0.0, ys = 0.0, xs = 0.0; + float sample_x = 0.0, sample_y = 0.0; + int x1 = 0, y1 = 0, sample_step = 0, pattern_size = 0; + int x2 = 0, y2 = 0, kx = 0, ky = 0, i = 0, j = 0, dcount = 0; + float fx = 0.0, fy = 0.0, ratio = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0, res4 = 0.0; + int scale = 0, dsize = 0, level = 0; + + // Subregion centers for the 4x4 gaussian weighting + float cx = -0.5, cy = 0.5; + + // Set the descriptor size and the sample and pattern sizes + dsize = 64; + sample_step = 5; + pattern_size = 12; + + // Get the information from the keypoint + ratio = (float)(1<(y1)+x1); + res2 = *(evolution_[level].Lx.ptr(y1)+x2); + res3 = *(evolution_[level].Lx.ptr(y2)+x1); + res4 = *(evolution_[level].Lx.ptr(y2)+x2); + rx = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Ly.ptr(y1)+x1); + res2 = *(evolution_[level].Ly.ptr(y1)+x2); + res3 = *(evolution_[level].Ly.ptr(y2)+x1); + res4 = *(evolution_[level].Ly.ptr(y2)+x2); + ry = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + rx = gauss_s1*rx; + ry = gauss_s1*ry; + + // Sum the derivatives to the cumulative descriptor + dx += rx; + dy += ry; + mdx += fabs(rx); + mdy += fabs(ry); + } + } + + // Add the values to the descriptor vector + gauss_s2 = gaussian(cx-2.0f,cy-2.0f,1.5f); + + desc[dcount++] = dx*gauss_s2; + desc[dcount++] = dy*gauss_s2; + desc[dcount++] = mdx*gauss_s2; + desc[dcount++] = mdy*gauss_s2; + + len += (dx*dx + dy*dy + mdx*mdx + mdy*mdy)*gauss_s2*gauss_s2; + + j += 9; + } + + i += 9; + } + + // convert to unit vector + len = sqrt(len); + + for (int i = 0; i < dsize; i++) { + desc[i] /= len; + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the descriptor of the provided keypoint given the + * main orientation of the keypoint + * @param kpt Input keypoint + * @param desc Descriptor vector + * @note Rectangular grid of 24 s x 24 s. Descriptor Length 64. The descriptor is inspired + * from Agrawal et al., CenSurE: Center Surround Extremas for Realtime Feature Detection and Matching, + * ECCV 2008 +*/ +void AKAZE::Get_MSURF_Descriptor_64(const cv::KeyPoint& kpt, float *desc) { + + float dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0, gauss_s1 = 0.0, gauss_s2 = 0.0; + 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; + float sample_x = 0.0, sample_y = 0.0, co = 0.0, si = 0.0, angle = 0.0; + float fx = 0.0, fy = 0.0, ratio = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0, res4 = 0.0; + int x1 = 0, y1 = 0, x2 = 0, y2 = 0, sample_step = 0, pattern_size = 0; + int kx = 0, ky = 0, i = 0, j = 0, dcount = 0; + int scale = 0, dsize = 0, level = 0; + + // Subregion centers for the 4x4 gaussian weighting + float cx = -0.5, cy = 0.5; + + // Set the descriptor size and the sample and pattern sizes + dsize = 64; + sample_step = 5; + pattern_size = 12; + + // Get the information from the keypoint + ratio = (float)(1<(y1)+x1); + res2 = *(evolution_[level].Lx.ptr(y1)+x2); + res3 = *(evolution_[level].Lx.ptr(y2)+x1); + res4 = *(evolution_[level].Lx.ptr(y2)+x2); + rx = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Ly.ptr(y1)+x1); + res2 = *(evolution_[level].Ly.ptr(y1)+x2); + res3 = *(evolution_[level].Ly.ptr(y2)+x1); + res4 = *(evolution_[level].Ly.ptr(y2)+x2); + ry = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + // Get the x and y derivatives on the rotated axis + rry = gauss_s1*(rx*co + ry*si); + rrx = gauss_s1*(-rx*si + ry*co); + + // Sum the derivatives to the cumulative descriptor + dx += rrx; + dy += rry; + mdx += fabs(rrx); + mdy += fabs(rry); + } + } + + // Add the values to the descriptor vector + gauss_s2 = gaussian(cx-2.0f,cy-2.0f,1.5f); + desc[dcount++] = dx*gauss_s2; + desc[dcount++] = dy*gauss_s2; + desc[dcount++] = mdx*gauss_s2; + desc[dcount++] = mdy*gauss_s2; + + len += (dx*dx + dy*dy + mdx*mdx + mdy*mdy)*gauss_s2*gauss_s2; + + j += 9; + } + + i += 9; + } + + // convert to unit vector + len = sqrt(len); + + for (int i = 0; i < dsize; i++) { + desc[i] /= len; + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the rupright descriptor (not rotation invariant) of + * the provided keypoint + * @param kpt Input keypoint + * @param desc Descriptor vector +*/ +void AKAZE::Get_Upright_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned char *desc) { + + float di = 0.0, dx = 0.0, dy = 0.0; + float ri = 0.0, rx = 0.0, ry = 0.0, xf = 0.0, yf = 0.0; + float sample_x = 0.0, sample_y = 0.0, ratio = 0.0; + int x1 = 0, y1 = 0, sample_step = 0, pattern_size = 0; + int level = 0, nsamples = 0, scale = 0; + int dcount1 = 0, dcount2 = 0; + + // Matrices for the M-LDB descriptor + Mat values_1 = Mat::zeros(4,descriptor_channels_,CV_32FC1); + Mat values_2 = Mat::zeros(9,descriptor_channels_,CV_32FC1); + Mat values_3 = Mat::zeros(16,descriptor_channels_,CV_32FC1); + + // Get the information from the keypoint + ratio = (float)(1<(y1)+x1); + rx = *(evolution_[level].Lx.ptr(y1)+x1); + ry = *(evolution_[level].Ly.ptr(y1)+x1); + + di += ri; + dx += rx; + dy += ry; + nsamples++; + } + } + + di /= nsamples; + dx /= nsamples; + dy /= nsamples; + + *(values_1.ptr(dcount2)) = di; + *(values_1.ptr(dcount2)+1) = dx; + *(values_1.ptr(dcount2)+2) = dy; + dcount2++; + } + } + + // Do binary comparison first level + for(int i = 0; i < 4; i++) { + for (int j = i+1; j < 4; j++) { + if (*(values_1.ptr(i)) > *(values_1.ptr(j))) { + desc[dcount1/8] |= (1<<(dcount1%8)); + } + dcount1++; + + if (*(values_1.ptr(i)+1) > *(values_1.ptr(j)+1)) { + desc[dcount1/8] |= (1<<(dcount1%8)); + } + dcount1++; + + if (*(values_1.ptr(i)+2) > *(values_1.ptr(j)+2)) { + desc[dcount1/8] |= (1<<(dcount1%8)); + } + dcount1++; + } + } + + // Second 3x3 grid + sample_step = ceil(pattern_size*2./3.); + dcount2 = 0; + + for (int i = -pattern_size; i < pattern_size; i+=sample_step) { + for (int j = -pattern_size; j < pattern_size; j+=sample_step) { + di=dx=dy=0.0; + nsamples = 0; + + for (int k = i; k < i + sample_step; k++) { + for (int l = j; l < j + sample_step; l++) { + // Get the coordinates of the sample point + sample_y = yf + l*scale; + sample_x = xf + k*scale; + + y1 = fRound(sample_y); + x1 = fRound(sample_x); + + ri = *(evolution_[level].Lt.ptr(y1)+x1); + rx = *(evolution_[level].Lx.ptr(y1)+x1); + ry = *(evolution_[level].Ly.ptr(y1)+x1); + + di += ri; + dx += rx; + dy += ry; + nsamples++; + } + } + + di /= nsamples; + dx /= nsamples; + dy /= nsamples; + + *(values_2.ptr(dcount2)) = di; + *(values_2.ptr(dcount2)+1) = dx; + *(values_2.ptr(dcount2)+2) = dy; + dcount2++; + } + } + + dcount2 = 0; + //Do binary comparison second level + for (int i = 0; i < 9; i++) { + for (int j = i+1; j < 9; j++) { + if (*(values_2.ptr(i)) > *(values_2.ptr(j))) { + desc[dcount1/8] |= (1<<(dcount1%8)); + } + dcount1++; + + if (*(values_2.ptr(i)+1) > *(values_2.ptr(j)+1)) { + desc[dcount1/8] |= (1<<(dcount1%8)); + } + dcount1++; + + if(*(values_2.ptr(i)+2) > *(values_2.ptr(j)+2)) { + desc[dcount1/8] |= (1<<(dcount1%8)); + } + dcount1++; + } + } + + // Third 4x4 grid + sample_step = pattern_size/2; + dcount2 = 0; + + for (int i = -pattern_size; i < pattern_size; i+=sample_step) { + for (int j = -pattern_size; j < pattern_size; j+=sample_step) { + di=dx=dy=0.0; + nsamples = 0; + + for (int k = i; k < i + sample_step; k++) { + for (int l = j; l < j + sample_step; l++) { + // Get the coordinates of the sample point + sample_y = yf + l*scale; + sample_x = xf + k*scale; + + y1 = fRound(sample_y); + x1 = fRound(sample_x); + + ri = *(evolution_[level].Lt.ptr(y1)+x1); + rx = *(evolution_[level].Lx.ptr(y1)+x1); + ry = *(evolution_[level].Ly.ptr(y1)+x1); + + di += ri; + dx += rx; + dy += ry; + nsamples++; + } + } + + di /= nsamples; + dx /= nsamples; + dy /= nsamples; + + *(values_3.ptr(dcount2)) = di; + *(values_3.ptr(dcount2)+1) = dx; + *(values_3.ptr(dcount2)+2) = dy; + dcount2++; + } + } + + dcount2 = 0; + //Do binary comparison third level + for (int i = 0; i < 16; i++) { + for (int j = i+1; j < 16; j++) { + if (*(values_3.ptr(i)) > *(values_3.ptr(j))) { + desc[dcount1/8] |= (1<<(dcount1%8)); + } + dcount1++; + + if (*(values_3.ptr(i)+1) > *(values_3.ptr(j)+1)) { + desc[dcount1/8] |= (1<<(dcount1%8)); + } + dcount1++; + + if (*(values_3.ptr(i)+2) > *(values_3.ptr(j)+2)) { + desc[dcount1/8] |= (1<<(dcount1%8)); + } + dcount1++; + } + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the descriptor of the provided keypoint given the + * main orientation of the keypoint + * @param kpt Input keypoint + * @param desc Descriptor vector +*/ +void AKAZE::Get_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned char *desc) { + + float di = 0.0, dx = 0.0, dy = 0.0, ratio = 0.0; + float ri = 0.0, rx = 0.0, ry = 0.0, rrx = 0.0, rry = 0.0, xf = 0.0, yf = 0.0; + float sample_x = 0.0, sample_y = 0.0, co = 0.0, si = 0.0, angle = 0.0; + int x1 = 0, y1 = 0, sample_step = 0, pattern_size = 0; + int level = 0, nsamples = 0, scale = 0; + int dcount1 = 0, dcount2 = 0; + + // Matrices for the M-LDB descriptor + Mat values_1 = Mat::zeros(4,descriptor_channels_,CV_32FC1); + Mat values_2 = Mat::zeros(9,descriptor_channels_,CV_32FC1); + Mat values_3 = Mat::zeros(16,descriptor_channels_,CV_32FC1); + + // Get the information from the keypoint + ratio = (float)(1<(y1)+x1); + rx = *(evolution_[level].Lx.ptr(y1)+x1); + ry = *(evolution_[level].Ly.ptr(y1)+x1); + + di += ri; + + if (descriptor_channels_ == 2) { + dx += sqrtf(rx*rx + ry*ry); + } + else if (descriptor_channels_ == 3) { + // Get the x and y derivatives on the rotated axis + rry = rx*co + ry*si; + rrx = -rx*si + ry*co; + dx += rrx; + dy += rry; + } + + nsamples++; + } + } + + di /= nsamples; + dx /= nsamples; + dy /= nsamples; + + *(values_1.ptr(dcount2)) = di; + if ( descriptor_channels_ > 1 ) { + *(values_1.ptr(dcount2)+1) = dx; + } + + if ( descriptor_channels_ > 2 ) { + *(values_1.ptr(dcount2)+2) = dy; + } + + dcount2++; + } + } + + // Do binary comparison first level + for (int i = 0; i < 4; i++) { + for (int j = i+1; j < 4; j++) { + if (*(values_1.ptr(i)) > *(values_1.ptr(j))) { + desc[dcount1/8] |= (1<<(dcount1%8)); + } + dcount1++; + } + } + + if (descriptor_channels_ > 1) { + for (int i = 0; i < 4; i++) { + for (int j = i+1; j < 4; j++) { + if (*(values_1.ptr(i)+1) > *(values_1.ptr(j)+1)) { + desc[dcount1/8] |= (1<<(dcount1%8)); + } + + dcount1++; + } + } + } + + if (descriptor_channels_ > 2) { + for (int i = 0; i < 4; i++) { + for ( int j = i+1; j < 4; j++) { + if (*(values_1.ptr(i)+2) > *(values_1.ptr(j)+2)) { + desc[dcount1/8] |= (1<<(dcount1%8)); + } + dcount1++; + } + } + } + + // Second 3x3 grid + sample_step = ceil(pattern_size*2./3.); + dcount2 = 0; + + for (int i = -pattern_size; i < pattern_size; i+=sample_step) { + for (int j = -pattern_size; j < pattern_size; j+=sample_step) { + + di=dx=dy=0.0; + nsamples = 0; + + for (int k = i; k < i + sample_step; k++) { + for (int l = j; l < j + sample_step; l++) { + + // Get the coordinates of the sample point + sample_y = yf + (l*scale*co + k*scale*si); + sample_x = xf + (-l*scale*si + k*scale*co); + + y1 = fRound(sample_y); + x1 = fRound(sample_x); + + ri = *(evolution_[level].Lt.ptr(y1)+x1); + rx = *(evolution_[level].Lx.ptr(y1)+x1); + ry = *(evolution_[level].Ly.ptr(y1)+x1); + di += ri; + + if (descriptor_channels_ == 2) { + dx += sqrtf(rx*rx + ry*ry); + } + else if (descriptor_channels_ == 3) { + // Get the x and y derivatives on the rotated axis + rry = rx*co + ry*si; + rrx = -rx*si + ry*co; + dx += rrx; + dy += rry; + } + + nsamples++; + } + } + + di /= nsamples; + dx /= nsamples; + dy /= nsamples; + + *(values_2.ptr(dcount2)) = di; + if (descriptor_channels_ > 1) { + *(values_2.ptr(dcount2)+1) = dx; + } + + if (descriptor_channels_ > 2) { + *(values_2.ptr(dcount2)+2) = dy; + } + + dcount2++; + } + } + + // Do binary comparison second level + for (int i = 0; i < 9; i++) { + for (int j = i+1; j < 9; j++) { + if (*(values_2.ptr(i)) > *(values_2.ptr(j))) { + desc[dcount1/8] |= (1<<(dcount1%8)); + } + dcount1++; + } + } + + if (descriptor_channels_ > 1) { + for (int i = 0; i < 9; i++) { + for (int j = i+1; j < 9; j++) { + if (*(values_2.ptr(i)+1) > *(values_2.ptr(j)+1)) { + desc[dcount1/8] |= (1<<(dcount1%8)); + } + dcount1++; + } + } + } + + if (descriptor_channels_ > 2) { + for (int i = 0; i < 9; i++) { + for (int j = i+1; j < 9; j++) { + if (*(values_2.ptr(i)+2) > *(values_2.ptr(j)+2)) { + desc[dcount1/8] |= (1<<(dcount1%8)); + } + dcount1++; + } + } + } + + // Third 4x4 grid + sample_step = pattern_size/2; + dcount2 = 0; + + for (int i = -pattern_size; i < pattern_size; i+=sample_step) { + for (int j = -pattern_size; j < pattern_size; j+=sample_step) { + di=dx=dy=0.0; + nsamples = 0; + + for (int k = i; k < i + sample_step; k++) { + for (int l = j; l < j + sample_step; l++) { + // Get the coordinates of the sample point + sample_y = yf + (l*scale*co + k*scale*si); + sample_x = xf + (-l*scale*si + k*scale*co); + + y1 = fRound(sample_y); + x1 = fRound(sample_x); + + ri = *(evolution_[level].Lt.ptr(y1)+x1); + rx = *(evolution_[level].Lx.ptr(y1)+x1); + ry = *(evolution_[level].Ly.ptr(y1)+x1); + di += ri; + + if (descriptor_channels_ == 2) { + dx += sqrtf(rx*rx + ry*ry); + } + else if (descriptor_channels_ == 3) { + // Get the x and y derivatives on the rotated axis + rry = rx*co + ry*si; + rrx = -rx*si + ry*co; + dx += rrx; + dy += rry; + } + + nsamples++; + } + } + + di /= nsamples; + dx /= nsamples; + dy /= nsamples; + + *(values_3.ptr(dcount2)) = di; + if (descriptor_channels_ > 1) { + *(values_3.ptr(dcount2)+1) = dx; + } + + if (descriptor_channels_ > 2) { + *(values_3.ptr(dcount2)+2) = dy; + } + + dcount2++; + } + } + + // Do binary comparison third level + for(int i = 0; i < 16; i++) { + for(int j = i+1; j < 16; j++) { + if (*(values_3.ptr(i)) > *(values_3.ptr(j))) { + desc[dcount1/8] |= (1<<(dcount1%8)); + } + dcount1++; + } + } + + if (descriptor_channels_ > 1) { + for (int i = 0; i < 16; i++) { + for (int j = i+1; j < 16; j++) { + if (*(values_3.ptr(i)+1) > *(values_3.ptr(j)+1)) { + desc[dcount1/8] |= (1<<(dcount1%8)); + } + dcount1++; + } + } + } + + if (descriptor_channels_ > 2) + { + for (int i = 0; i < 16; i++) { + for (int j = i+1; j < 16; j++) { + if (*(values_3.ptr(i)+2) > *(values_3.ptr(j)+2)) { + desc[dcount1/8] |= (1<<(dcount1%8)); + } + dcount1++; + } + } + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the M-LDB descriptor of the provided keypoint given the + * main orientation of the keypoint. The descriptor is computed based on a subset of + * the bits of the whole descriptor + * @param kpt Input keypoint + * @param desc Descriptor vector +*/ +void AKAZE::Get_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char *desc) { + + float di, dx, dy; + float rx, ry; + float sample_x = 0.f, sample_y = 0.f; + int x1 = 0, y1 = 0; + + // Get the information from the keypoint + float ratio = (float)(1<::zeros((4+9+16)*descriptor_channels_,1); + + // Sample everything, but only do the comparisons + vector steps(3); + steps.at(0) = descriptor_pattern_size_; + steps.at(1) = ceil(2.f*descriptor_pattern_size_/3.f); + steps.at(2) = descriptor_pattern_size_/2; + + for (int i=0; i(i); + int sample_step = steps.at(coords[0]); + di=0.0f; + dx=0.0f; + dy=0.0f; + + for (int k = coords[1]; k < coords[1] + sample_step; k++) { + for (int l = coords[2]; l < coords[2] + sample_step; l++) { + // Get the coordinates of the sample point + sample_y = yf + (l*scale*co + k*scale*si); + sample_x = xf + (-l*scale*si + k*scale*co); + + y1 = fRound(sample_y); + x1 = fRound(sample_x); + + di += *(evolution_[level].Lt.ptr(y1)+x1); + + if (descriptor_channels_ > 1) { + rx = *(evolution_[level].Lx.ptr(y1)+x1); + ry = *(evolution_[level].Ly.ptr(y1)+x1); + + if (descriptor_channels_ == 2) { + dx += sqrtf(rx*rx + ry*ry); + } + else if (descriptor_channels_ == 3) { + // Get the x and y derivatives on the rotated axis + dx += rx*co + ry*si; + dy += -rx*si + ry*co; + } + } + } + } + + *(values.ptr(descriptor_channels_*i)) = di; + + if (descriptor_channels_ == 2) { + *(values.ptr(descriptor_channels_*i+1)) = dx; + } + else if (descriptor_channels_ == 3) { + *(values.ptr(descriptor_channels_*i+1)) = dx; + *(values.ptr(descriptor_channels_*i+2)) = dy; + } + } + + // Do the comparisons + const float *vals = values.ptr(0); + const int *comps = descriptorBits_.ptr(0); + + for (int i=0; i vals[comps[2*i +1]]) { + desc[i/8] |= (1<<(i%8)); + } + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the upright (not rotation invariant) M-LDB descriptor + * of the provided keypoint given the main orientation of the keypoint. + * The descriptor is computed based on a subset of the bits of the whole descriptor + * @param kpt Input keypoint + * @param desc Descriptor vector +*/ +void AKAZE::Get_Upright_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char *desc) { + + float di = 0.0f, dx = 0.0f, dy = 0.0f; + float rx = 0.0f, ry = 0.0f; + float sample_x = 0.0f, sample_y = 0.0f; + int x1 = 0, y1 = 0; + + // Get the information from the keypoint + float ratio = (float)(1<::zeros((4+9+16)*descriptor_channels_,1); + + vector steps(3); + steps.at(0) = descriptor_pattern_size_; + steps.at(1) = ceil(2.f*descriptor_pattern_size_/3.f); + steps.at(2) = descriptor_pattern_size_/2; + + for (int i=0; i < descriptorSamples_.rows; i++) { + int *coords = descriptorSamples_.ptr(i); + int sample_step = steps.at(coords[0]); + di=0.0f; + dx=0.0f; + dy=0.0f; + + for (int k = coords[1]; k < coords[1] + sample_step; k++) { + for (int l = coords[2]; l < coords[2] + sample_step; l++) { + // Get the coordinates of the sample point + sample_y = yf + l*scale; + sample_x = xf + k*scale; + + y1 = fRound(sample_y); + x1 = fRound(sample_x); + di += *(evolution_[level].Lt.ptr(y1)+x1); + + if (descriptor_channels_ > 1) { + rx = *(evolution_[level].Lx.ptr(y1)+x1); + ry = *(evolution_[level].Ly.ptr(y1)+x1); + + if (descriptor_channels_ == 2) { + dx += sqrtf(rx*rx + ry*ry); + } + else if (descriptor_channels_ == 3) { + dx += rx; + dy += ry; + } + } + } + } + + *(values.ptr(descriptor_channels_*i)) = di; + + if (descriptor_channels_ == 2) { + *(values.ptr(descriptor_channels_*i+1)) = dx; + } + else if (descriptor_channels_ == 3) { + *(values.ptr(descriptor_channels_*i+1)) = dx; + *(values.ptr(descriptor_channels_*i+2)) = dy; + } + } + + // Do the comparisons + const float *vals = values.ptr(0); + const int *comps = descriptorBits_.ptr(0); + + for (int i=0; i vals[comps[2*i +1]]) { + desc[i/8] |= (1<<(i%8)); + } + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method displays the computation times +*/ +void AKAZE::Show_Computation_Times(void) { + + cout << "(*) Time Scale Space: " << tscale_ << endl; + cout << "(*) Time Detector: " << tdetector_ << endl; + cout << " - Time Derivatives: " << tderivatives_ << endl; + cout << " - Time Extrema: " << textrema_ << endl; + cout << " - Time Subpixel: " << tsubpixel_ << endl; + cout << "(*) Time Descriptor: " << tdescriptor_ << endl; +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function computes a (quasi-random) list of bits to be taken + * from the full descriptor. To speed the extraction, the function creates + * a list of the samples that are involved in generating at least a bit (sampleList) + * and a list of the comparisons between those samples (comparisons) + * @param sampleList + * @param comparisons The matrix with the binary comparisons + * @param nbits The number of bits of the descriptor + * @param pattern_size The pattern size for the binary descriptor + * @param nchannels Number of channels to consider in the descriptor (1-3) + * @note The function keeps the 18 bits (3-channels by 6 comparisons) of the + * coarser grid, since it provides the most robust estimations + */ +void generateDescriptorSubsample(cv::Mat& sampleList, cv::Mat& comparisons, int nbits, + int pattern_size, int nchannels) { + + int ssz = 0; + for (int i=0; i<3; i++) { + int gz = (i+2)*(i+2); + ssz += gz*(gz-1)/2; + } + ssz *= nchannels; + + CV_Assert(nbits<=ssz && "descriptor size can't be bigger than full descriptor"); + + // Since the full descriptor is usually under 10k elements, we pick + // the selection from the full matrix. We take as many samples per + // pick as the number of channels. For every pick, we + // take the two samples involved and put them in the sampling list + + Mat_ fullM(ssz/nchannels,5); + for (size_t i=0, c=0; i<3; i++) { + int gdiv = i+2; //grid divisions, per row + int gsz = gdiv*gdiv; + int psz = ceil(2.*pattern_size/(float)gdiv); + + for (int j=0; j comps = Mat_(nchannels*ceil(nbits/(float)nchannels),2); + comps = 1000; + + // Select some samples. A sample includes all channels + int count =0; + size_t npicks = ceil(nbits/(float)nchannels); + Mat_ samples(29,3); + Mat_ fullcopy = fullM.clone(); + samples = -1; + + for (size_t i=0; i= 0 && y >= 0) { + return atanf(y/x); + } + + if (x < 0 && y >= 0) { + return CV_PI - atanf(-y/x); + } + + if (x < 0 && y < 0) { + return CV_PI + atanf(y/x); + } + + if(x >= 0 && y < 0) { + return 2.0*CV_PI - atanf(-y/x); + } + + return 0; +} + +//************************************************************************************** +//************************************************************************************** + +/** + * @brief This function computes the value of a 2D Gaussian function + * @param x X Position + * @param y Y Position + * @param sig Standard Deviation +*/ +inline float gaussian(float x, float y, float sigma) { + + return expf(-(x*x+y*y)/(2.0f*sigma*sigma)); +} + +//************************************************************************************** +//************************************************************************************** + +/** + * @brief This function checks descriptor limits + * @param x X Position + * @param y Y Position + * @param width Image width + * @param height Image height +*/ +inline void check_descriptor_limits(int &x, int &y, const int width, const int height) { + + if (x < 0) { + x = 0; + } + + if (y < 0) { + y = 0; + } + + if (x > width-1) { + x = width-1; + } + + if (y > height-1) { + y = height-1; + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This funtion rounds float to nearest integer + * @param flt Input float + * @return dst Nearest integer + */ +inline int fRound(float flt) +{ + return (int)(flt+0.5f); +} + diff --git a/modules/features2d/src/akaze/AKAZE.h b/modules/features2d/src/akaze/AKAZE.h new file mode 100644 index 0000000000..fd1ec07fa3 --- /dev/null +++ b/modules/features2d/src/akaze/AKAZE.h @@ -0,0 +1,175 @@ +/** + * @file AKAZE.h + * @brief Main class for detecting and computing binary descriptors in an + * accelerated nonlinear scale space + * @date Mar 27, 2013 + * @author Pablo F. Alcantarilla, Jesus Nuevo + */ + +#ifndef _AKAZE_H_ +#define _AKAZE_H_ + +//************************************************************************************* +//************************************************************************************* + +// Includes +#include "config.h" +#include "fed.h" +#include "utils.h" +#include "nldiffusion_functions.h" + +//************************************************************************************* +//************************************************************************************* + +// AKAZE Class Declaration +class AKAZE { + +private: + + // Parameters of the AKAZE class + int omax_; // Maximum octave level + int noctaves_; // Number of octaves + int nsublevels_; // Number of sublevels per octave level + int img_width_; // Width of the original image + int img_height_; // Height of the original image + float soffset_; // Base scale offset + float factor_size_; // Factor for the multiscale derivatives + float sderivatives_; // Standard deviation of the Gaussian for the nonlinear diff. derivatives + float kcontrast_; // The contrast parameter for the scalar nonlinear diffusion + float dthreshold_; // Feature detector threshold response + int diffusivity_; // Diffusivity type, 0->PM G1, 1->PM G2, 2-> Weickert, 3->Charbonnier + int descriptor_; // Descriptor mode: + // 0-> SURF_UPRIGHT, 1->SURF + // 2-> M-SURF_UPRIGHT, 3->M-SURF + // 4-> M-LDB_UPRIGHT, 5->M-LDB + int descriptor_size_; // Size of the descriptor in bits. Use 0 for the full descriptor + int descriptor_pattern_size_; // Size of the pattern. Actual size sampled is 2*pattern_size + int descriptor_channels_; // Number of channels to consider in the M-LDB descriptor + bool save_scale_space_; // For saving scale space images + bool verbosity_; // Verbosity level + std::vector evolution_; // Vector of nonlinear diffusion evolution + + // FED parameters + int ncycles_; // Number of cycles + bool reordering_; // Flag for reordering time steps + std::vector > tsteps_; // Vector of FED dynamic time steps + std::vector nsteps_; // Vector of number of steps per cycle + + // Some matrices for the M-LDB descriptor computation + cv::Mat descriptorSamples_; // List of positions in the grids to sample LDB bits from. + cv::Mat descriptorBits_; + cv::Mat bitMask_; + + // Computation times variables in ms + double tkcontrast_; // Kcontrast factor computation + double tscale_; // Nonlinear Scale space generation + double tderivatives_; // Multiscale derivatives + double tdetector_; // Feature detector + double textrema_; // Scale Space extrema + double tsubpixel_; // Subpixel refinement + double tdescriptor_; // Feature descriptors + +public: + + // Constructor + AKAZE(const AKAZEOptions &options); + + // Destructor + ~AKAZE(void); + + // Setters + void Set_Octave_Max(const int& omax) { + omax_ = omax; + } + void Set_NSublevels(const int& nsublevels) { + nsublevels_ = nsublevels; + } + void Set_Save_Scale_Space_Flag(const bool& save_scale_space) { + save_scale_space_ = save_scale_space; + } + void Set_Image_Width(const int& img_width) { + img_width_ = img_width; + } + void Set_Image_Height(const int& img_height) { + img_height_ = img_height; + } + + // Getters + int Get_Image_Width(void) { + return img_width_; + } + int Get_Image_Height(void) { + return img_height_; + } + double Get_Time_KContrast(void) { + return tkcontrast_; + } + double Get_Time_Scale_Space(void) { + return tscale_; + } + double Get_Time_Derivatives(void) { + return tderivatives_; + } + double Get_Time_Detector(void) { + return tdetector_; + } + double Get_Time_Descriptor(void) { + return tdescriptor_; + } + + // Scale Space methods + void Allocate_Memory_Evolution(void); + int Create_Nonlinear_Scale_Space(const cv::Mat& img); + void Feature_Detection(std::vector& kpts); + void Compute_Determinant_Hessian_Response(void); + void Compute_Multiscale_Derivatives(void); + void Find_Scale_Space_Extrema(std::vector& kpts); + void Do_Subpixel_Refinement(std::vector& kpts); + void Feature_Suppression_Distance(std::vector& kpts, float mdist); + + // Feature description methods + void Compute_Descriptors(std::vector& kpts, cv::Mat& desc); + void Compute_Main_Orientation_SURF(cv::KeyPoint& kpt); + + // SURF Pattern Descriptor + void Get_SURF_Descriptor_Upright_64(const cv::KeyPoint& kpt, float *desc); + void Get_SURF_Descriptor_64(const cv::KeyPoint& kpt, float *desc); + + // M-SURF Pattern Descriptor + void Get_MSURF_Upright_Descriptor_64(const cv::KeyPoint& kpt, float *desc); + void Get_MSURF_Descriptor_64(const cv::KeyPoint& kpt, float *desc); + + // M-LDB Pattern Descriptor + void Get_Upright_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned char *desc); + void Get_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned char *desc); + void Get_Upright_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char *desc); + void Get_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char *desc); + + // Methods for saving some results and showing computation times + void Save_Scale_Space(void); + void Save_Detector_Responses(void); + void Show_Computation_Times(void); +}; + +//************************************************************************************* +//************************************************************************************* + +// Inline functions +/** + * @brief This function sets default parameters for the A-KAZE detector. + * @param options AKAZE options + */ +void setDefaultAKAZEOptions(AKAZEOptions& options); + +// Inline functions +void generateDescriptorSubsample(cv::Mat& sampleList, cv::Mat& comparisons, + int nbits, int pattern_size, int nchannels); +float get_angle(float x, float y); +float gaussian(float x, float y, float sigma); +void check_descriptor_limits(int& x, int& y, const int width, const int height); +int fRound(float flt); + +//************************************************************************************* +//************************************************************************************* + +#endif diff --git a/modules/features2d/src/akaze/config.h b/modules/features2d/src/akaze/config.h new file mode 100644 index 0000000000..331c89275f --- /dev/null +++ b/modules/features2d/src/akaze/config.h @@ -0,0 +1,155 @@ +#ifndef _CONFIG_H_ +#define _CONFIG_H_ + +// STL +#include +#include +#include +#include +#include + +// OpenCV +#include "precomp.hpp" + +// OpenMP +#ifdef _OPENMP +# include +#endif + +// Lookup table for 2d gaussian (sigma = 2.5) where (0,0) is top left and (6,6) is bottom right +const float gauss25[7][7] = { + {0.02546481f, 0.02350698f, 0.01849125f, 0.01239505f, 0.00708017f, 0.00344629f, 0.00142946f}, + {0.02350698f, 0.02169968f, 0.01706957f, 0.01144208f, 0.00653582f, 0.00318132f, 0.00131956f}, + {0.01849125f, 0.01706957f, 0.01342740f, 0.00900066f, 0.00514126f, 0.00250252f, 0.00103800f}, + {0.01239505f, 0.01144208f, 0.00900066f, 0.00603332f, 0.00344629f, 0.00167749f, 0.00069579f}, + {0.00708017f, 0.00653582f, 0.00514126f, 0.00344629f, 0.00196855f, 0.00095820f, 0.00039744f}, + {0.00344629f, 0.00318132f, 0.00250252f, 0.00167749f, 0.00095820f, 0.00046640f, 0.00019346f}, + {0.00142946f, 0.00131956f, 0.00103800f, 0.00069579f, 0.00039744f, 0.00019346f, 0.00008024f} +}; + + +// Scale Space parameters +const float DEFAULT_SCALE_OFFSET = 1.60f; // Base scale offset (sigma units) +const float DEFAULT_FACTOR_SIZE = 1.5f; // Factor for the multiscale derivatives +const int DEFAULT_OCTAVE_MIN = 0; // Initial octave level (-1 means that the size of the input image is duplicated) +const int DEFAULT_OCTAVE_MAX = 4; // Maximum octave evolution of the image 2^sigma (coarsest scale sigma units) +const int DEFAULT_NSUBLEVELS = 4; // Default number of sublevels per scale level +const int DEFAULT_DIFFUSIVITY_TYPE = 1; +const float KCONTRAST_PERCENTILE = 0.7f; +const int KCONTRAST_NBINS = 300; +const float DEFAULT_SIGMA_SMOOTHING_DERIVATIVES = 1.0f; +const float DEFAULT_KCONTRAST = .01f; + + +// Detector Parameters +const float DEFAULT_DETECTOR_THRESHOLD = 0.001f; // Detector response threshold to accept point +const float DEFAULT_MIN_DETECTOR_THRESHOLD = 0.00001f; // Minimum Detector response threshold to accept point +const int DEFAULT_LDB_DESCRIPTOR_SIZE = 0; // Use 0 for the full descriptor, or the number of bits +const int DEFAULT_LDB_PATTERN_SIZE = 10; // Actual patch size is 2*pattern_size*point.scale; +const int DEFAULT_LDB_CHANNELS = 3; + +// Descriptor Parameters +enum DESCRIPTOR_TYPE +{ + SURF_UPRIGHT = 0, // Upright descriptors, not invariant to rotation + SURF = 1, + MSURF_UPRIGHT = 2, // Upright descriptors, not invariant to rotation + MSURF = 3, + MLDB_UPRIGHT = 4, // Upright descriptors, not invariant to rotation + MLDB = 5 +}; + +const int DEFAULT_DESCRIPTOR = MLDB; + +// Some debugging options +const bool DEFAULT_SAVE_SCALE_SPACE = false; // For saving the scale space images +const bool DEFAULT_VERBOSITY = false; // Verbosity level (0->no verbosity) +const bool DEFAULT_SHOW_RESULTS = true; // For showing the output image with the detected features plus some ratios +const bool DEFAULT_SAVE_KEYPOINTS = false; // For saving the list of keypoints + +// Options structure +struct AKAZEOptions +{ + int omin; + int omax; + int nsublevels; + int img_width; + int img_height; + int diffusivity; + float soffset; + float sderivatives; + float dthreshold; + float dthreshold2; + int descriptor; + int descriptor_size; + int descriptor_channels; + int descriptor_pattern_size; + bool save_scale_space; + bool save_keypoints; + bool verbosity; + + AKAZEOptions() + { + // Load the default options + soffset = DEFAULT_SCALE_OFFSET; + omax = DEFAULT_OCTAVE_MAX; + nsublevels = DEFAULT_NSUBLEVELS; + dthreshold = DEFAULT_DETECTOR_THRESHOLD; + diffusivity = DEFAULT_DIFFUSIVITY_TYPE; + descriptor = DEFAULT_DESCRIPTOR; + descriptor_size = DEFAULT_LDB_DESCRIPTOR_SIZE; + descriptor_channels = DEFAULT_LDB_CHANNELS; + descriptor_pattern_size = DEFAULT_LDB_PATTERN_SIZE; + sderivatives = DEFAULT_SIGMA_SMOOTHING_DERIVATIVES; + save_scale_space = DEFAULT_SAVE_SCALE_SPACE; + save_keypoints = DEFAULT_SAVE_KEYPOINTS; + verbosity = DEFAULT_VERBOSITY; + } + + friend std::ostream& operator<<(std::ostream& os, + const AKAZEOptions& akaze_options) + { + os << std::left; +#define CHECK_AKAZE_OPTION(option) \ + os << std::setw(33) << #option << " = " << option << std::endl + + // Scale-space parameters. + CHECK_AKAZE_OPTION(akaze_options.omax); + CHECK_AKAZE_OPTION(akaze_options.nsublevels); + CHECK_AKAZE_OPTION(akaze_options.soffset); + CHECK_AKAZE_OPTION(akaze_options.sderivatives); + CHECK_AKAZE_OPTION(akaze_options.diffusivity); + // Detection parameters. + CHECK_AKAZE_OPTION(akaze_options.dthreshold); + // Descriptor parameters. + CHECK_AKAZE_OPTION(akaze_options.descriptor); + CHECK_AKAZE_OPTION(akaze_options.descriptor_channels); + CHECK_AKAZE_OPTION(akaze_options.descriptor_size); + // Save scale-space + CHECK_AKAZE_OPTION(akaze_options.save_scale_space); + // Verbose option for debug. + CHECK_AKAZE_OPTION(akaze_options.verbosity); +#undef CHECK_AKAZE_OPTIONS + + return os; + } +}; + +struct tevolution +{ + cv::Mat Lx, Ly; // First order spatial derivatives + cv::Mat Lxx, Lxy, Lyy; // Second order spatial derivatives + cv::Mat Lflow; // Diffusivity image + cv::Mat Lt; // Evolution image + cv::Mat Lsmooth; // Smoothed image + cv::Mat Lstep; // Evolution step update + cv::Mat Ldet; // Detector response + float etime; // Evolution time + float esigma; // Evolution sigma. For linear diffusion t = sigma^2 / 2 + int octave; // Image octave + int sublevel; // Image sublevel in each octave + int sigma_size; // Integer sigma. For computing the feature detector responses +}; + + +#endif \ No newline at end of file diff --git a/modules/features2d/src/akaze/fed.h b/modules/features2d/src/akaze/fed.h new file mode 100644 index 0000000000..4ac82f68e3 --- /dev/null +++ b/modules/features2d/src/akaze/fed.h @@ -0,0 +1,26 @@ +#ifndef FED_H +#define FED_H + +//****************************************************************************** +//****************************************************************************** + +// Includes +#include +#include + +//************************************************************************************* +//************************************************************************************* + +// Declaration of functions +int fed_tau_by_process_time(const float& T, const int& M, const float& tau_max, + const bool& reordering, std::vector& tau); +int fed_tau_by_cycle_time(const float& t, const float& tau_max, + const bool& reordering, std::vector &tau) ; +int fed_tau_internal(const int& n, const float& scale, const float& tau_max, + const bool& reordering, std::vector &tau); +bool fed_is_prime_internal(const int& number); + +//************************************************************************************* +//************************************************************************************* + +#endif // FED_H diff --git a/modules/features2d/src/akaze/nldiffusion_functions.cpp b/modules/features2d/src/akaze/nldiffusion_functions.cpp new file mode 100644 index 0000000000..0699e92ca0 --- /dev/null +++ b/modules/features2d/src/akaze/nldiffusion_functions.cpp @@ -0,0 +1,431 @@ +//============================================================================= +// +// nldiffusion_functions.cpp +// Authors: Pablo F. Alcantarilla (1), Jesus Nuevo (2) +// Institutions: Georgia Institute of Technology (1) +// TrueVision Solutions (2) +// Date: 15/09/2013 +// Email: pablofdezalc@gmail.com +// +// AKAZE Features Copyright 2013, Pablo F. Alcantarilla, Jesus Nuevo +// All Rights Reserved +// See LICENSE for the license information +//============================================================================= + +/** + * @file nldiffusion_functions.cpp + * @brief Functions for nonlinear diffusion filtering applications + * @date Sep 15, 2013 + * @author Pablo F. Alcantarilla, Jesus Nuevo + */ + +#include "nldiffusion_functions.h" + +using namespace std; +using namespace cv; + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function smoothes an image with a Gaussian kernel + * @param src Input image + * @param dst Output image + * @param ksize_x Kernel size in X-direction (horizontal) + * @param ksize_y Kernel size in Y-direction (vertical) + * @param sigma Kernel standard deviation + */ +void gaussian_2D_convolution(const cv::Mat& src, cv::Mat& dst, const size_t& ksize_x, + const size_t& ksize_y, const float& sigma) { + + size_t ksize_x_ = 0, ksize_y_ = 0; + + // Compute an appropriate kernel size according to the specified sigma + if (sigma > ksize_x || sigma > ksize_y || ksize_x == 0 || ksize_y == 0) { + ksize_x_ = ceil(2.0*(1.0 + (sigma-0.8)/(0.3))); + ksize_y_ = ksize_x_; + } + + // The kernel size must be and odd number + if ((ksize_x_ % 2) == 0) { + ksize_x_ += 1; + } + + if ((ksize_y_ % 2) == 0) { + ksize_y_ += 1; + } + + // Perform the Gaussian Smoothing with border replication + GaussianBlur(src,dst,Size(ksize_x_,ksize_y_),sigma,sigma,BORDER_REPLICATE); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function computes image derivatives with Scharr kernel + * @param src Input image + * @param dst Output image + * @param xorder Derivative order in X-direction (horizontal) + * @param yorder Derivative order in Y-direction (vertical) + * @note Scharr operator approximates better rotation invariance than + * other stencils such as Sobel. See Weickert and Scharr, + * A Scheme for Coherence-Enhancing Diffusion Filtering with Optimized Rotation Invariance, + * Journal of Visual Communication and Image Representation 2002 + */ +void image_derivatives_scharr(const cv::Mat& src, cv::Mat& dst, + const size_t& xorder, const size_t& yorder) { + Scharr(src,dst,CV_32F,xorder,yorder,1.0,0,BORDER_DEFAULT); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function computes the Perona and Malik conductivity coefficient g1 + * g1 = exp(-|dL|^2/k^2) + * @param Lx First order image derivative in X-direction (horizontal) + * @param Ly First order image derivative in Y-direction (vertical) + * @param dst Output image + * @param k Contrast factor parameter + */ +void pm_g1(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k) { + exp(-(Lx.mul(Lx)+Ly.mul(Ly))/(k*k),dst); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function computes the Perona and Malik conductivity coefficient g2 + * g2 = 1 / (1 + dL^2 / k^2) + * @param Lx First order image derivative in X-direction (horizontal) + * @param Ly First order image derivative in Y-direction (vertical) + * @param dst Output image + * @param k Contrast factor parameter + */ +void pm_g2(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k) { + dst = 1.0/(1.0+(Lx.mul(Lx)+Ly.mul(Ly))/(k*k)); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function computes Weickert conductivity coefficient gw + * @param Lx First order image derivative in X-direction (horizontal) + * @param Ly First order image derivative in Y-direction (vertical) + * @param dst Output image + * @param k Contrast factor parameter + * @note For more information check the following paper: J. Weickert + * Applications of nonlinear diffusion in image processing and computer vision, + * Proceedings of Algorithmy 2000 + */ +void weickert_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k) { + Mat modg; + pow((Lx.mul(Lx) + Ly.mul(Ly))/(k*k),4,modg); + cv::exp(-3.315/modg, dst); + dst = 1.0 - dst; +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function computes Charbonnier conductivity coefficient gc + * gc = 1 / sqrt(1 + dL^2 / k^2) + * @param Lx First order image derivative in X-direction (horizontal) + * @param Ly First order image derivative in Y-direction (vertical) + * @param dst Output image + * @param k Contrast factor parameter + * @note For more information check the following paper: J. Weickert + * Applications of nonlinear diffusion in image processing and computer vision, + * Proceedings of Algorithmy 2000 + */ +void charbonnier_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k) { + Mat den; + cv::sqrt(1.0+(Lx.mul(Lx)+Ly.mul(Ly))/(k*k),den); + dst = 1.0/ den; +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function computes a good empirical value for the k contrast factor + * given an input image, the percentile (0-1), the gradient scale and the number of + * bins in the histogram + * @param img Input image + * @param perc Percentile of the image gradient histogram (0-1) + * @param gscale Scale for computing the image gradient histogram + * @param nbins Number of histogram bins + * @param ksize_x Kernel size in X-direction (horizontal) for the Gaussian smoothing kernel + * @param ksize_y Kernel size in Y-direction (vertical) for the Gaussian smoothing kernel + * @return k contrast factor + */ +float compute_k_percentile(const cv::Mat& img, const float& perc, const float& gscale, + const size_t& nbins, const size_t& ksize_x, const size_t& ksize_y) { + + size_t nbin = 0, nelements = 0, nthreshold = 0, k = 0; + float kperc = 0.0, modg = 0.0, lx = 0.0, ly = 0.0; + float npoints = 0.0; + float hmax = 0.0; + + // Create the array for the histogram + float *hist = new float[nbins]; + + // Create the matrices + Mat gaussian = Mat::zeros(img.rows,img.cols,CV_32F); + Mat Lx = Mat::zeros(img.rows,img.cols,CV_32F); + Mat Ly = Mat::zeros(img.rows,img.cols,CV_32F); + + // Set the histogram to zero, just in case + for (size_t i = 0; i < nbins; i++) { + hist[i] = 0.0; + } + + // Perform the Gaussian convolution + gaussian_2D_convolution(img,gaussian,ksize_x,ksize_y,gscale); + + // Compute the Gaussian derivatives Lx and Ly + image_derivatives_scharr(gaussian,Lx,1,0); + image_derivatives_scharr(gaussian,Ly,0,1); + + // Skip the borders for computing the histogram + for (int i = 1; i < gaussian.rows-1; i++) { + for (int j = 1; j < gaussian.cols-1; j++) { + lx = *(Lx.ptr(i)+j); + ly = *(Ly.ptr(i)+j); + modg = sqrt(lx*lx + ly*ly); + + // Get the maximum + if (modg > hmax) { + hmax = modg; + } + } + } + + // Skip the borders for computing the histogram + for (int i = 1; i < gaussian.rows-1; i++) { + for (int j = 1; j < gaussian.cols-1; j++) { + lx = *(Lx.ptr(i)+j); + ly = *(Ly.ptr(i)+j); + modg = sqrt(lx*lx + ly*ly); + + // Find the correspondent bin + if (modg != 0.0) { + nbin = floor(nbins*(modg/hmax)); + + if (nbin == nbins) { + nbin--; + } + + hist[nbin]++; + npoints++; + } + } + } + + // Now find the perc of the histogram percentile + nthreshold = (size_t)(npoints*perc); + + for (k = 0; nelements < nthreshold && k < nbins; k++) { + nelements = nelements + hist[k]; + } + + if (nelements < nthreshold) { + kperc = 0.03; + } + else { + kperc = hmax*((float)(k)/(float)nbins); + } + + delete [] hist; + return kperc; +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function computes Scharr image derivatives + * @param src Input image + * @param dst Output image + * @param xorder Derivative order in X-direction (horizontal) + * @param yorder Derivative order in Y-direction (vertical) + * @param scale Scale factor for the derivative size + */ +void compute_scharr_derivatives(const cv::Mat& src, cv::Mat& dst, const size_t& xorder, + const size_t& yorder, const size_t& scale) { + + Mat kx, ky; + compute_derivative_kernels(kx, ky, xorder,yorder,scale); + sepFilter2D(src,dst,CV_32F,kx,ky); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function performs a scalar non-linear diffusion step + * @param Ld2 Output image in the evolution + * @param c Conductivity image + * @param Lstep Previous image in the evolution + * @param stepsize The step size in time units + * @note Forward Euler Scheme 3x3 stencil + * The function c is a scalar value that depends on the gradient norm + * dL_by_ds = d(c dL_by_dx)_by_dx + d(c dL_by_dy)_by_dy + */ +void nld_step_scalar(cv::Mat& Ld, const cv::Mat& c, cv::Mat& Lstep, const float& stepsize) { + +#ifdef _OPENMP +#pragma omp parallel for schedule(dynamic) +#endif + for (int i = 1; i < Lstep.rows-1; i++) { + for (int j = 1; j < Lstep.cols-1; j++) { + float xpos = ((*(c.ptr(i)+j))+(*(c.ptr(i)+j+1)))*((*(Ld.ptr(i)+j+1))-(*(Ld.ptr(i)+j))); + float xneg = ((*(c.ptr(i)+j-1))+(*(c.ptr(i)+j)))*((*(Ld.ptr(i)+j))-(*(Ld.ptr(i)+j-1))); + + float ypos = ((*(c.ptr(i)+j))+(*(c.ptr(i+1)+j)))*((*(Ld.ptr(i+1)+j))-(*(Ld.ptr(i)+j))); + float yneg = ((*(c.ptr(i-1)+j))+(*(c.ptr(i)+j)))*((*(Ld.ptr(i)+j))-(*(Ld.ptr(i-1)+j))); + + *(Lstep.ptr(i)+j) = 0.5*stepsize*(xpos-xneg + ypos-yneg); + } + } + + for (int j = 1; j < Lstep.cols-1; j++) { + float xpos = ((*(c.ptr(0)+j))+(*(c.ptr(0)+j+1)))*((*(Ld.ptr(0)+j+1))-(*(Ld.ptr(0)+j))); + float xneg = ((*(c.ptr(0)+j-1))+(*(c.ptr(0)+j)))*((*(Ld.ptr(0)+j))-(*(Ld.ptr(0)+j-1))); + + float ypos = ((*(c.ptr(0)+j))+(*(c.ptr(1)+j)))*((*(Ld.ptr(1)+j))-(*(Ld.ptr(0)+j))); + float yneg = ((*(c.ptr(0)+j))+(*(c.ptr(0)+j)))*((*(Ld.ptr(0)+j))-(*(Ld.ptr(0)+j))); + + *(Lstep.ptr(0)+j) = 0.5*stepsize*(xpos-xneg + ypos-yneg); + } + + for (int j = 1; j < Lstep.cols-1; j++) { + float xpos = ((*(c.ptr(Lstep.rows-1)+j))+(*(c.ptr(Lstep.rows-1)+j+1)))*((*(Ld.ptr(Lstep.rows-1)+j+1))-(*(Ld.ptr(Lstep.rows-1)+j))); + float xneg = ((*(c.ptr(Lstep.rows-1)+j-1))+(*(c.ptr(Lstep.rows-1)+j)))*((*(Ld.ptr(Lstep.rows-1)+j))-(*(Ld.ptr(Lstep.rows-1)+j-1))); + + float ypos = ((*(c.ptr(Lstep.rows-1)+j))+(*(c.ptr(Lstep.rows-1)+j)))*((*(Ld.ptr(Lstep.rows-1)+j))-(*(Ld.ptr(Lstep.rows-1)+j))); + float yneg = ((*(c.ptr(Lstep.rows-2)+j))+(*(c.ptr(Lstep.rows-1)+j)))*((*(Ld.ptr(Lstep.rows-1)+j))-(*(Ld.ptr(Lstep.rows-2)+j))); + + *(Lstep.ptr(Lstep.rows-1)+j) = 0.5*stepsize*(xpos-xneg + ypos-yneg); + } + + for (int i = 1; i < Lstep.rows-1; i++) { + float xpos = ((*(c.ptr(i)))+(*(c.ptr(i)+1)))*((*(Ld.ptr(i)+1))-(*(Ld.ptr(i)))); + float xneg = ((*(c.ptr(i)))+(*(c.ptr(i))))*((*(Ld.ptr(i)))-(*(Ld.ptr(i)))); + + float ypos = ((*(c.ptr(i)))+(*(c.ptr(i+1))))*((*(Ld.ptr(i+1)))-(*(Ld.ptr(i)))); + float yneg = ((*(c.ptr(i-1)))+(*(c.ptr(i))))*((*(Ld.ptr(i)))-(*(Ld.ptr(i-1)))); + + *(Lstep.ptr(i)) = 0.5*stepsize*(xpos-xneg + ypos-yneg); + } + + for (int i = 1; i < Lstep.rows-1; i++) { + float xpos = ((*(c.ptr(i)+Lstep.cols-1))+(*(c.ptr(i)+Lstep.cols-1)))*((*(Ld.ptr(i)+Lstep.cols-1))-(*(Ld.ptr(i)+Lstep.cols-1))); + float xneg = ((*(c.ptr(i)+Lstep.cols-2))+(*(c.ptr(i)+Lstep.cols-1)))*((*(Ld.ptr(i)+Lstep.cols-1))-(*(Ld.ptr(i)+Lstep.cols-2))); + + float ypos = ((*(c.ptr(i)+Lstep.cols-1))+(*(c.ptr(i+1)+Lstep.cols-1)))*((*(Ld.ptr(i+1)+Lstep.cols-1))-(*(Ld.ptr(i)+Lstep.cols-1))); + float yneg = ((*(c.ptr(i-1)+Lstep.cols-1))+(*(c.ptr(i)+Lstep.cols-1)))*((*(Ld.ptr(i)+Lstep.cols-1))-(*(Ld.ptr(i-1)+Lstep.cols-1))); + + *(Lstep.ptr(i)+Lstep.cols-1) = 0.5*stepsize*(xpos-xneg + ypos-yneg); + } + + Ld = Ld + Lstep; +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function downsamples the input image with the kernel [1/4,1/2,1/4] + * @param img Input image to be downsampled + * @param dst Output image with half of the resolution of the input image + */ +void downsample_image(const cv::Mat& src, cv::Mat& dst) { + + int i1 = 0, j1 = 0, i2 = 0, j2 = 0; + + for (i1 = 1; i1 < src.rows; i1+=2) { + j2 = 0; + for (j1 = 1; j1 < src.cols; j1+=2) { + *(dst.ptr(i2)+j2) = 0.5*(*(src.ptr(i1)+j1))+0.25*(*(src.ptr(i1)+j1-1) + *(src.ptr(i1)+j1+1)); + j2++; + } + + i2++; + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function downsamples the input image using OpenCV resize + * @param img Input image to be downsampled + * @param dst Output image with half of the resolution of the input image + */ +void halfsample_image(const cv::Mat& src, cv::Mat& dst) { + + // Make sure the destination image is of the right size + CV_Assert(src.cols/2==dst.cols); + CV_Assert(src.rows / 2 == dst.rows); + resize(src,dst,dst.size(),0,0,cv::INTER_AREA); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief Compute Scharr derivative kernels for sizes different than 3 + * @param kx_ The derivative kernel in x-direction + * @param ky_ The derivative kernel in y-direction + * @param dx The derivative order in x-direction + * @param dy The derivative order in y-direction + * @param scale The kernel size + */ +void compute_derivative_kernels(cv::OutputArray kx_, cv::OutputArray ky_, + const size_t& dx, const size_t& dy, const size_t& scale) { + + const int ksize = 3 + 2*(scale-1); + + // The usual Scharr kernel + if (scale == 1) { + getDerivKernels(kx_,ky_,dx,dy,0,true,CV_32F); + return; + } + + kx_.create(ksize,1,CV_32F,-1,true); + ky_.create(ksize,1,CV_32F,-1,true); + Mat kx = kx_.getMat(); + Mat ky = ky_.getMat(); + + float w = 10.0/3.0; + float norm = 1.0/(2.0*scale*(w+2.0)); + + for (int k = 0; k < 2; k++) { + Mat* kernel = k == 0 ? &kx : &ky; + int order = k == 0 ? dx : dy; + float kerI[1000]; + + for (int t = 0; trows, kernel->cols, CV_32F, &kerI[0]); + temp.copyTo(*kernel); + } +} diff --git a/modules/features2d/src/akaze/nldiffusion_functions.h b/modules/features2d/src/akaze/nldiffusion_functions.h new file mode 100644 index 0000000000..172fa25f3e --- /dev/null +++ b/modules/features2d/src/akaze/nldiffusion_functions.h @@ -0,0 +1,41 @@ +#ifndef _NLDIFFUSION_FUNCTIONS_H_ +#define _NLDIFFUSION_FUNCTIONS_H_ + +//****************************************************************************** +//****************************************************************************** + +// Includes +#include "precomp.hpp" + +// OpenMP Includes +#ifdef _OPENMP +# include +#endif + +//************************************************************************************* +//************************************************************************************* + +// Declaration of functions +void gaussian_2D_convolution(const cv::Mat& src, cv::Mat& dst, const size_t& ksize_x, + const size_t& ksize_y, const float& sigma); +void image_derivatives_scharr(const cv::Mat& src, cv::Mat& dst, + const size_t& xorder, const size_t& yorder); +void pm_g1(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k); +void pm_g2(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k); +void weickert_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k); +void charbonnier_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k); +float compute_k_percentile(const cv::Mat& img, const float& perc, const float& gscale, + const size_t& nbins, const size_t& ksize_x, const size_t& ksize_y); +void compute_scharr_derivatives(const cv::Mat& src, cv::Mat& dst, const size_t& xorder, + const size_t& yorder, const size_t& scale); +void nld_step_scalar(cv::Mat& Ld, const cv::Mat& c, cv::Mat& Lstep, const float& stepsize); +void downsample_image(const cv::Mat& src, cv::Mat& dst); +void halfsample_image(const cv::Mat& src, cv::Mat& dst); +void compute_derivative_kernels(cv::OutputArray kx_, cv::OutputArray ky_, + const size_t& dx, const size_t& dy, const size_t& scale); + +//************************************************************************************* +//************************************************************************************* + + +#endif diff --git a/modules/features2d/src/akaze/utils.cpp b/modules/features2d/src/akaze/utils.cpp new file mode 100644 index 0000000000..eb14abcd51 --- /dev/null +++ b/modules/features2d/src/akaze/utils.cpp @@ -0,0 +1,196 @@ +//============================================================================= +// +// utils.cpp +// Authors: Pablo F. Alcantarilla (1), Jesus Nuevo (2) +// Institutions: Georgia Institute of Technology (1) +// TrueVision Solutions (2) +// +// Date: 15/09/2013 +// Email: pablofdezalc@gmail.com +// +// AKAZE Features Copyright 2013, Pablo F. Alcantarilla, Jesus Nuevo +// All Rights Reserved +// See LICENSE for the license information +//============================================================================= + +/** + * @file utils.cpp + * @brief Some utilities functions + * @date Sep 15, 2013 + * @author Pablo F. Alcantarilla, Jesus Nuevo + */ + +#include "precomp.hpp" +#include "utils.h" + +// Namespaces +using namespace std; +using namespace cv; + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function computes the minimum value of a float image + * @param src Input image + * @param value Minimum value + */ +void compute_min_32F(const cv::Mat &src, float &value) { + + float aux = 1000.0; + + for (int i = 0; i < src.rows; i++) { + for (int j = 0; j < src.cols; j++) { + if (src.at(i,j) < aux) { + aux = src.at(i,j); + } + } + } + + value = aux; +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function computes the maximum value of a float image + * @param src Input image + * @param value Maximum value + */ +void compute_max_32F(const cv::Mat &src, float &value) { + + float aux = 0.0; + + for (int i = 0; i < src.rows; i++) { + for (int j = 0; j < src.cols; j++) { + if (src.at(i,j) > aux) { + aux = src.at(i,j); + } + } + } + + value = aux; +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function converts the scale of the input image prior to visualization + * @param src Input/Output image + * @param value Maximum value + */ +void convert_scale(cv::Mat &src) { + + float min_val = 0, max_val = 0; + + compute_min_32F(src,min_val); + + src = src - min_val; + + compute_max_32F(src,max_val); + src = src / max_val; +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function copies the input image and converts the scale of the copied + * image prior visualization + * @param src Input image + * @param dst Output image + */ +void copy_and_convert_scale(const cv::Mat &src, cv::Mat dst) { + + float min_val = 0, max_val = 0; + + src.copyTo(dst); + compute_min_32F(dst,min_val); + + dst = dst - min_val; + + compute_max_32F(dst,max_val); + dst = dst / max_val; +} + +//************************************************************************************* +//************************************************************************************* + +const size_t length = string("--descriptor_channels").size() + 2; +static inline std::ostream& cout_help() +{ cout << setw(length); return cout; } + +static inline std::string toUpper(std::string s) +{ + std::transform(s.begin(), s.end(), s.begin(), ::toupper); + return s; +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function shows the possible command line configuration options + */ +void show_input_options_help(int example) { + + fflush(stdout); + cout << "A-KAZE Features" << endl; + cout << "Usage: "; + if (example == 0) { + cout << "./akaze_features -i img.jpg [options]" << endl; + } + else if (example == 1) { + cout << "./akaze_match img1.jpg img2.pgm homography.txt [options]" << endl; + } + else if (example == 2) { + cout << "./akaze_compare img1.jpg img2.pgm homography.txt [options]" << endl; + } + + cout << endl; + cout_help() << "Options below are not mandatory. Unless specified, default arguments are used." << endl << endl; + // Justify on the left + cout << left; + // Generalities + cout_help() << "--help" << "Show the command line options" << endl; + cout_help() << "--verbose " << "Verbosity is required" << endl; + cout_help() << endl; + // Scale-space parameters + cout_help() << "--soffset" << "Base scale offset (sigma units)" << endl; + cout_help() << "--omax" << "Maximum octave of image evolution" << endl; + cout_help() << "--nsublevels" << "Number of sublevels per octave" << endl; + cout_help() << "--diffusivity" << "Diffusivity function. Possible values:" << endl; + cout_help() << " " << "0 -> Perona-Malik, g1 = exp(-|dL|^2/k^2)" << endl; + cout_help() << " " << "1 -> Perona-Malik, g2 = 1 / (1 + dL^2 / k^2)" << endl; + cout_help() << " " << "2 -> Weickert diffusivity" << endl; + cout_help() << " " << "3 -> Charbonnier diffusivity" << endl; + cout_help() << endl; + // Feature detection parameters. + cout_help() << "--dthreshold" << "Feature detector threshold response for keypoints" << endl; + cout_help() << " " << "(0.001 can be a good value)" << endl; + cout_help() << endl; + // Descriptor parameters. + cout_help() << "--descriptor" << "Descriptor Type. Possible values:" << endl; + cout_help() << " " << "0 -> SURF_UPRIGHT" << endl; + cout_help() << " " << "1 -> SURF" << endl; + cout_help() << " " << "2 -> M-SURF_UPRIGHT," << endl; + cout_help() << " " << "3 -> M-SURF" << endl; + cout_help() << " " << "4 -> M-LDB_UPRIGHT" << endl; + cout_help() << " " << "5 -> M-LDB" << endl; + + cout_help() << "--descriptor_channels " << "Descriptor Channels for M-LDB. Valid values: " << endl; + cout_help() << " " << "1 -> intensity" << endl; + cout_help() << " " << "2 -> intensity + gradient magnitude" << endl; + cout_help() << " " << "3 -> intensity + X and Y gradients" < show detection results." << endl; + cout_help() << " " << "0 -> don't show detection results" << endl; + cout_help() << endl; +} diff --git a/modules/features2d/src/akaze/utils.h b/modules/features2d/src/akaze/utils.h new file mode 100644 index 0000000000..894c836ed0 --- /dev/null +++ b/modules/features2d/src/akaze/utils.h @@ -0,0 +1,54 @@ + +#ifndef _UTILS_H_ +#define _UTILS_H_ + +//****************************************************************************** +//****************************************************************************** + +// OpenCV Includes +#include "precomp.hpp" + +// System Includes +#include +#include +#include +#include +#include +#include +#include + +//****************************************************************************** +//****************************************************************************** + +// Stringify common types such as int, double and others. +template +inline std::string to_string(const T& x) { + std::stringstream oss; + oss << x; + return oss.str(); +} + +//****************************************************************************** +//****************************************************************************** + +// Stringify and format integral types as follows: +// to_formatted_string( 1, 2) produces string: '01' +// to_formatted_string( 5, 2) produces string: '05' +// to_formatted_string( 19, 2) produces string: '19' +// to_formatted_string( 19, 3) produces string: '019' +template +inline std::string to_formatted_string(Integer x, int num_digits) { + std::stringstream oss; + oss << std::setfill('0') << std::setw(num_digits) << x; + return oss.str(); +} + +//****************************************************************************** +//****************************************************************************** + +void compute_min_32F(const cv::Mat& src, float& value); +void compute_max_32F(const cv::Mat& src, float& value); +void convert_scale(cv::Mat& src); +void copy_and_convert_scale(const cv::Mat& src, cv::Mat& dst); + +#endif diff --git a/modules/features2d/src/kaze.cpp b/modules/features2d/src/kaze.cpp new file mode 100644 index 0000000000..e69de29bb2 diff --git a/modules/features2d/src/kaze/KAZE.cpp b/modules/features2d/src/kaze/KAZE.cpp new file mode 100644 index 0000000000..09246e1670 --- /dev/null +++ b/modules/features2d/src/kaze/KAZE.cpp @@ -0,0 +1,2801 @@ + +//============================================================================= +// +// KAZE.cpp +// Author: Pablo F. Alcantarilla +// Institution: University d'Auvergne +// Address: Clermont Ferrand, France +// Date: 21/01/2012 +// Email: pablofdezalc@gmail.com +// +// KAZE Features Copyright 2012, Pablo F. Alcantarilla +// All Rights Reserved +// See LICENSE for the license information +//============================================================================= + +/** + * @file KAZE.cpp + * @brief Main class for detecting and describing features in a nonlinear + * scale space + * @date Jan 21, 2012 + * @author Pablo F. Alcantarilla + */ + +#include "KAZE.h" + +// Namespaces +using namespace std; +using namespace cv; + +//******************************************************************************* +//******************************************************************************* + +/** + * @brief KAZE constructor with input options + * @param options KAZE configuration options + * @note The constructor allocates memory for the nonlinear scale space +*/ +KAZE::KAZE(KAZEOptions& options) { + + soffset_ = options.soffset; + sderivatives_ = options.sderivatives; + omax_ = options.omax; + nsublevels_ = options.nsublevels; + save_scale_space_ = options.save_scale_space; + verbosity_ = options.verbosity; + img_width_ = options.img_width; + img_height_ = options.img_height; + dthreshold_ = options.dthreshold; + diffusivity_ = options.diffusivity; + descriptor_mode_ = options.descriptor; + use_fed_ = options.use_fed; + use_upright_ = options.upright; + use_extended_ = options.extended; + kcontrast_ = DEFAULT_KCONTRAST; + ncycles_ = 0; + reordering_ = true; + tkcontrast_ = 0.0; + tnlscale_ = 0.0; + tdetector_ = 0.0; + tmderivatives_ = 0.0; + tdresponse_ = 0.0; + tdescriptor_ = 0.0; + + // Now allocate memory for the evolution + Allocate_Memory_Evolution(); +} + +//******************************************************************************* +//******************************************************************************* + +/** + * @brief KAZE destructor +*/ +KAZE::~KAZE(void) { + + evolution_.clear(); +} + +//******************************************************************************* +//******************************************************************************* + +/** + * @brief This method allocates the memory for the nonlinear diffusion evolution +*/ +void KAZE::Allocate_Memory_Evolution(void) { + + // Allocate the dimension of the matrices for the evolution + for (int i = 0; i <= omax_-1; i++) { + for (int j = 0; j <= nsublevels_-1; j++) { + + TEvolution aux; + aux.Lx = cv::Mat::zeros(img_height_,img_width_,CV_32F); + aux.Ly = cv::Mat::zeros(img_height_,img_width_,CV_32F); + aux.Lxx = cv::Mat::zeros(img_height_,img_width_,CV_32F); + aux.Lxy = cv::Mat::zeros(img_height_,img_width_,CV_32F); + aux.Lyy = cv::Mat::zeros(img_height_,img_width_,CV_32F); + aux.Lflow = cv::Mat::zeros(img_height_,img_width_,CV_32F); + aux.Lt = cv::Mat::zeros(img_height_,img_width_,CV_32F); + aux.Lsmooth = cv::Mat::zeros(img_height_,img_width_,CV_32F); + aux.Lstep = cv::Mat::zeros(img_height_,img_width_,CV_32F); + aux.Ldet = cv::Mat::zeros(img_height_,img_width_,CV_32F); + aux.esigma = soffset_*pow((float)2.0,(float)(j)/(float)(nsublevels_) + i); + aux.etime = 0.5*(aux.esigma*aux.esigma); + aux.sigma_size = fRound(aux.esigma); + aux.octave = i; + aux.sublevel = j; + evolution_.push_back(aux); + } + } + + // Allocate memory for the FED number of cycles and time steps + if (use_fed_) { + for (size_t i = 1; i < evolution_.size(); i++) { + int naux = 0; + vector tau; + float ttime = 0.0; + ttime = evolution_[i].etime-evolution_[i-1].etime; + naux = fed_tau_by_process_time(ttime,1,0.25,reordering_,tau); + nsteps_.push_back(naux); + tsteps_.push_back(tau); + ncycles_++; + } + } + else { + // Allocate memory for the auxiliary variables that are used in the AOS scheme + Ltx_ = Mat::zeros(img_width_,img_height_,CV_32F); + Lty_ = Mat::zeros(img_height_,img_width_,CV_32F); + px_ = Mat::zeros(img_height_,img_width_,CV_32F); + py_ = Mat::zeros(img_height_,img_width_,CV_32F); + ax_ = Mat::zeros(img_height_,img_width_,CV_32F); + ay_ = Mat::zeros(img_height_,img_width_,CV_32F); + bx_ = Mat::zeros(img_height_-1,img_width_,CV_32F); + by_ = Mat::zeros(img_height_-1,img_width_,CV_32F); + qr_ = Mat::zeros(img_height_-1,img_width_,CV_32F); + qc_ = Mat::zeros(img_height_,img_width_-1,CV_32F); + } + +} + +//******************************************************************************* +//******************************************************************************* + +/** + * @brief This method creates the nonlinear scale space for a given image + * @param img Input image for which the nonlinear scale space needs to be created + * @return 0 if the nonlinear scale space was created successfully. -1 otherwise +*/ +int KAZE::Create_Nonlinear_Scale_Space(const cv::Mat &img) { + + double t2 = 0.0, t1 = 0.0; + + if (evolution_.size() == 0) { + cout << "Error generating the nonlinear scale space!!" << endl; + cout << "Firstly you need to call KAZE::Allocate_Memory_Evolution()" << endl; + return -1; + } + + t1 = getTickCount(); + + // Copy the original image to the first level of the evolution + img.copyTo(evolution_[0].Lt); + gaussian_2D_convolution(evolution_[0].Lt,evolution_[0].Lt,0,0,soffset_); + gaussian_2D_convolution(evolution_[0].Lt,evolution_[0].Lsmooth,0,0,sderivatives_); + + // Firstly compute the kcontrast factor + Compute_KContrast(evolution_[0].Lt,KCONTRAST_PERCENTILE); + + t2 = getTickCount(); + tkcontrast_ = 1000.0*(t2-t1) / getTickFrequency(); + + if (verbosity_ == true) { + cout << "Computed image evolution step. Evolution time: " << evolution_[0].etime << + " Sigma: " << evolution_[0].esigma << endl; + } + + // Now generate the rest of evolution levels + for ( size_t i = 1; i < evolution_.size(); i++) { + + evolution_[i-1].Lt.copyTo(evolution_[i].Lt); + gaussian_2D_convolution(evolution_[i-1].Lt,evolution_[i].Lsmooth,0,0,sderivatives_); + + // Compute the Gaussian derivatives Lx and Ly + Scharr(evolution_[i].Lsmooth,evolution_[i].Lx,CV_32F,1,0,1,0,BORDER_DEFAULT); + Scharr(evolution_[i].Lsmooth,evolution_[i].Ly,CV_32F,0,1,1,0,BORDER_DEFAULT); + + // Compute the conductivity equation + if (diffusivity_ == 0) { + pm_g1(evolution_[i].Lx,evolution_[i].Ly,evolution_[i].Lflow,kcontrast_); + } + else if (diffusivity_ == 1) { + pm_g2(evolution_[i].Lx,evolution_[i].Ly,evolution_[i].Lflow,kcontrast_); + } + else if (diffusivity_ == 2) { + weickert_diffusivity(evolution_[i].Lx,evolution_[i].Ly,evolution_[i].Lflow,kcontrast_); + } + + // Perform FED n inner steps + if (use_fed_) { + for (int j = 0; j < nsteps_[i-1]; j++) { + nld_step_scalar(evolution_[i].Lt,evolution_[i].Lflow,evolution_[i].Lstep,tsteps_[i-1][j]); + } + } + else { + // Perform the evolution step with AOS + AOS_Step_Scalar(evolution_[i].Lt,evolution_[i-1].Lt,evolution_[i].Lflow, + evolution_[i].etime-evolution_[i-1].etime); + } + + if (verbosity_ == true) { + cout << "Computed image evolution step " << i << " Evolution time: " << evolution_[i].etime << + " Sigma: " << evolution_[i].esigma << endl; + } + } + + t2 = getTickCount(); + tnlscale_ = 1000.0*(t2-t1) / getTickFrequency(); + + return 0; +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the k contrast factor + * @param img Input image + * @param kpercentile Percentile of the gradient histogram +*/ +void KAZE::Compute_KContrast(const cv::Mat &img, const float &kpercentile) { + + if (verbosity_ == true) { + cout << "Computing Kcontrast factor." << endl; + } + + if (COMPUTE_KCONTRAST == true) { + kcontrast_ = compute_k_percentile(img,kpercentile,sderivatives_,KCONTRAST_NBINS,0,0); + } + + if (verbosity_ == true) { + cout << "kcontrast = " << kcontrast_ << endl; + cout << endl << "Now computing the nonlinear scale space!!" << endl; + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the multiscale derivatives for the nonlinear scale space +*/ +void KAZE::Compute_Multiscale_Derivatives(void) +{ + double t2 = 0.0, t1 = 0.0; + t1 = getTickCount(); + +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (size_t i = 0; i < evolution_.size(); i++) { + + if (verbosity_ == true) { + cout << "Computing multiscale derivatives. Evolution time: " << evolution_[i].etime + << " Step (pixels): " << evolution_[i].sigma_size << endl; + } + + // Compute multiscale derivatives for the detector + compute_scharr_derivatives(evolution_[i].Lsmooth,evolution_[i].Lx,1,0,evolution_[i].sigma_size); + compute_scharr_derivatives(evolution_[i].Lsmooth,evolution_[i].Ly,0,1,evolution_[i].sigma_size); + compute_scharr_derivatives(evolution_[i].Lx,evolution_[i].Lxx,1,0,evolution_[i].sigma_size); + compute_scharr_derivatives(evolution_[i].Ly,evolution_[i].Lyy,0,1,evolution_[i].sigma_size); + compute_scharr_derivatives(evolution_[i].Lx,evolution_[i].Lxy,0,1,evolution_[i].sigma_size); + + evolution_[i].Lx = evolution_[i].Lx*((evolution_[i].sigma_size)); + evolution_[i].Ly = evolution_[i].Ly*((evolution_[i].sigma_size)); + evolution_[i].Lxx = evolution_[i].Lxx*((evolution_[i].sigma_size)*(evolution_[i].sigma_size)); + evolution_[i].Lxy = evolution_[i].Lxy*((evolution_[i].sigma_size)*(evolution_[i].sigma_size)); + evolution_[i].Lyy = evolution_[i].Lyy*((evolution_[i].sigma_size)*(evolution_[i].sigma_size)); + } + + t2 = getTickCount(); + tmderivatives_ = 1000.0*(t2-t1) / getTickFrequency(); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the feature detector response for the nonlinear scale space + * @note We use the Hessian determinant as feature detector +*/ +void KAZE::Compute_Detector_Response(void) { + + double t2 = 0.0, t1 = 0.0; + float lxx = 0.0, lxy = 0.0, lyy = 0.0; + + t1 = getTickCount(); + + // Firstly compute the multiscale derivatives + Compute_Multiscale_Derivatives(); + + for (size_t i = 0; i < evolution_.size(); i++) { + + // Determinant of the Hessian + if (verbosity_ == true) { + cout << "Computing detector response. Determinant of Hessian. Evolution time: " << evolution_[i].etime << endl; + } + + for (int ix = 0; ix < img_height_; ix++) { + for (int jx = 0; jx < img_width_; jx++) { + lxx = *(evolution_[i].Lxx.ptr(ix)+jx); + lxy = *(evolution_[i].Lxy.ptr(ix)+jx); + lyy = *(evolution_[i].Lyy.ptr(ix)+jx); + *(evolution_[i].Ldet.ptr(ix)+jx) = (lxx*lyy-lxy*lxy); + } + } + } + + t2 = getTickCount(); + tdresponse_ = 1000.0*(t2-t1) / getTickFrequency(); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method selects interesting keypoints through the nonlinear scale space + * @param kpts Vector of keypoints +*/ +void KAZE::Feature_Detection(std::vector& kpts) { + + double t2 = 0.0, t1 = 0.0; + t1 = getTickCount(); + + // Firstly compute the detector response for each pixel and scale level + Compute_Detector_Response(); + + // Find scale space extrema + Determinant_Hessian_Parallel(kpts); + + // Perform some subpixel refinement + Do_Subpixel_Refinement(kpts); + + t2 = getTickCount(); + tdetector_ = 1000.0*(t2-t1) / getTickFrequency(); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method performs the detection of keypoints by using the normalized + * score of the Hessian determinant through the nonlinear scale space + * @param kpts Vector of keypoints + * @note We compute features for each of the nonlinear scale space level in a different processing thread +*/ +void KAZE::Determinant_Hessian_Parallel(std::vector& kpts) { + + int level = 0; + float dist = 0.0, smax = 3.0; + int npoints = 0, id_repeated = 0; + int left_x = 0, right_x = 0, up_y = 0, down_y = 0; + bool is_extremum = false, is_repeated = false, is_out = false; + + // Delete the memory of the vector of keypoints vectors + // In case we use the same kaze object for multiple images + for (size_t i = 0; i < kpts_par_.size(); i++) { + vector().swap(kpts_par_[i]); + } + kpts_par_.clear(); + vector aux; + + // Allocate memory for the vector of vectors + for (size_t i = 1; i < evolution_.size()-1; i++) { + kpts_par_.push_back(aux); + } + +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (size_t i = 1; i < evolution_.size()-1; i++) { + Find_Extremum_Threading(i); + } + + // Now fill the vector of keypoints!!! + for (size_t i = 0; i < kpts_par_.size(); i++) { + for (size_t j = 0; j < kpts_par_[i].size(); j++) { + level = i+1; + is_extremum = true; + is_repeated = false; + is_out = false; + + // Check in case we have the same point as maxima in previous evolution levels + for (size_t ik = 0; ik < kpts.size(); ik++) { + if (kpts[ik].class_id == level || kpts[ik].class_id == level+1 || kpts[ik].class_id == level-1) { + dist = pow(kpts_par_[i][j].pt.x-kpts[ik].pt.x,2)+pow(kpts_par_[i][j].pt.y-kpts[ik].pt.y,2); + + if (dist < evolution_[level].sigma_size*evolution_[level].sigma_size) { + if (kpts_par_[i][j].response > kpts[ik].response) { + id_repeated = ik; + is_repeated = true; + } + else { + is_extremum = false; + } + + break; + } + } + } + + if (is_extremum == true) { + // Check that the point is under the image limits for the descriptor computation + left_x = fRound(kpts_par_[i][j].pt.x-smax*kpts_par_[i][j].size); + right_x = fRound(kpts_par_[i][j].pt.x+smax*kpts_par_[i][j].size); + up_y = fRound(kpts_par_[i][j].pt.y-smax*kpts_par_[i][j].size); + down_y = fRound(kpts_par_[i][j].pt.y+smax*kpts_par_[i][j].size); + + if (left_x < 0 || right_x >= evolution_[level].Ldet.cols || + up_y < 0 || down_y >= evolution_[level].Ldet.rows) { + is_out = true; + } + + is_out = false; + + if (is_out == false) { + if (is_repeated == false) { + kpts.push_back(kpts_par_[i][j]); + npoints++; + } + else { + kpts[id_repeated] = kpts_par_[i][j]; + } + } + } + } + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method is called by the thread which is responsible of finding extrema + * at a given nonlinear scale level + * @param level Index in the nonlinear scale space evolution +*/ +void KAZE::Find_Extremum_Threading(const int& level) { + + float value = 0.0; + bool is_extremum = false; + + for (int ix = 1; ix < img_height_-1; ix++) { + for (int jx = 1; jx < img_width_-1; jx++) { + + is_extremum = false; + value = *(evolution_[level].Ldet.ptr(ix)+jx); + + // Filter the points with the detector threshold + if (value > dthreshold_ && value >= DEFAULT_MIN_DETECTOR_THRESHOLD) { + if (value >= *(evolution_[level].Ldet.ptr(ix)+jx-1)) { + // First check on the same scale + if (check_maximum_neighbourhood(evolution_[level].Ldet,1,value,ix,jx,1)) { + // Now check on the lower scale + if (check_maximum_neighbourhood(evolution_[level-1].Ldet,1,value,ix,jx,0)) { + // Now check on the upper scale + if (check_maximum_neighbourhood(evolution_[level+1].Ldet,1,value,ix,jx,0)) { + is_extremum = true; + } + } + } + } + } + + // Add the point of interest!! + if (is_extremum == true) { + KeyPoint point; + point.pt.x = jx; + point.pt.y = ix; + point.response = fabs(value); + point.size = evolution_[level].esigma; + point.octave = evolution_[level].octave; + point.class_id = level; + + // We use the angle field for the sublevel value + // Then, we will replace this angle field with the main orientation + point.angle = evolution_[level].sublevel; + kpts_par_[level-1].push_back(point); + } + } + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method performs subpixel refinement of the detected keypoints + * @param kpts Vector of detected keypoints +*/ +void KAZE::Do_Subpixel_Refinement(std::vector &kpts) { + + int step = 1; + int x = 0, y = 0; + float Dx = 0.0, Dy = 0.0, Ds = 0.0, dsc = 0.0; + float Dxx = 0.0, Dyy = 0.0, Dss = 0.0, Dxy = 0.0, Dxs = 0.0, Dys = 0.0; + Mat A = Mat::zeros(3,3,CV_32F); + Mat b = Mat::zeros(3,1,CV_32F); + Mat dst = Mat::zeros(3,1,CV_32F); + double t2 = 0.0, t1 = 0.0; + + t1 = cv::getTickCount(); + vector kpts_(kpts); + + for (size_t i = 0; i < kpts_.size(); i++) { + + x = kpts_[i].pt.x; + y = kpts_[i].pt.y; + + // Compute the gradient + Dx = (1.0/(2.0*step))*(*(evolution_[kpts_[i].class_id].Ldet.ptr(y)+x+step) + -*(evolution_[kpts_[i].class_id].Ldet.ptr(y)+x-step)); + Dy = (1.0/(2.0*step))*(*(evolution_[kpts_[i].class_id].Ldet.ptr(y+step)+x) + -*(evolution_[kpts_[i].class_id].Ldet.ptr(y-step)+x)); + Ds = 0.5*(*(evolution_[kpts_[i].class_id+1].Ldet.ptr(y)+x) + -*(evolution_[kpts_[i].class_id-1].Ldet.ptr(y)+x)); + + // Compute the Hessian + Dxx = (1.0/(step*step))*(*(evolution_[kpts_[i].class_id].Ldet.ptr(y)+x+step) + + *(evolution_[kpts_[i].class_id].Ldet.ptr(y)+x-step) + -2.0*(*(evolution_[kpts_[i].class_id].Ldet.ptr(y)+x))); + + Dyy = (1.0/(step*step))*(*(evolution_[kpts_[i].class_id].Ldet.ptr(y+step)+x) + + *(evolution_[kpts_[i].class_id].Ldet.ptr(y-step)+x) + -2.0*(*(evolution_[kpts_[i].class_id].Ldet.ptr(y)+x))); + + Dss = *(evolution_[kpts_[i].class_id+1].Ldet.ptr(y)+x) + + *(evolution_[kpts_[i].class_id-1].Ldet.ptr(y)+x) + -2.0*(*(evolution_[kpts_[i].class_id].Ldet.ptr(y)+x)); + + Dxy = (1.0/(4.0*step))*(*(evolution_[kpts_[i].class_id].Ldet.ptr(y+step)+x+step) + +(*(evolution_[kpts_[i].class_id].Ldet.ptr(y-step)+x-step))) + -(1.0/(4.0*step))*(*(evolution_[kpts_[i].class_id].Ldet.ptr(y-step)+x+step) + +(*(evolution_[kpts_[i].class_id].Ldet.ptr(y+step)+x-step))); + + Dxs = (1.0/(4.0*step))*(*(evolution_[kpts_[i].class_id+1].Ldet.ptr(y)+x+step) + +(*(evolution_[kpts_[i].class_id-1].Ldet.ptr(y)+x-step))) + -(1.0/(4.0*step))*(*(evolution_[kpts_[i].class_id+1].Ldet.ptr(y)+x-step) + +(*(evolution_[kpts_[i].class_id-1].Ldet.ptr(y)+x+step))); + + Dys = (1.0/(4.0*step))*(*(evolution_[kpts_[i].class_id+1].Ldet.ptr(y+step)+x) + +(*(evolution_[kpts_[i].class_id-1].Ldet.ptr(y-step)+x))) + -(1.0/(4.0*step))*(*(evolution_[kpts_[i].class_id+1].Ldet.ptr(y-step)+x) + +(*(evolution_[kpts_[i].class_id-1].Ldet.ptr(y+step)+x))); + + // Solve the linear system + *(A.ptr(0)) = Dxx; + *(A.ptr(1)+1) = Dyy; + *(A.ptr(2)+2) = Dss; + + *(A.ptr(0)+1) = *(A.ptr(1)) = Dxy; + *(A.ptr(0)+2) = *(A.ptr(2)) = Dxs; + *(A.ptr(1)+2) = *(A.ptr(2)+1) = Dys; + + *(b.ptr(0)) = -Dx; + *(b.ptr(1)) = -Dy; + *(b.ptr(2)) = -Ds; + + solve(A,b,dst,DECOMP_LU); + + if (fabs(*(dst.ptr(0))) <= 1.0 && fabs(*(dst.ptr(1))) <= 1.0 && fabs(*(dst.ptr(2))) <= 1.0) { + kpts_[i].pt.x += *(dst.ptr(0)); + kpts_[i].pt.y += *(dst.ptr(1)); + dsc = kpts_[i].octave + (kpts_[i].angle+*(dst.ptr(2)))/((float)(nsublevels_)); + + // In OpenCV the size of a keypoint is the diameter!! + kpts_[i].size = 2.0*soffset_*pow((float)2.0,dsc); + kpts_[i].angle = 0.0; + } + // Set the points to be deleted after the for loop + else { + kpts_[i].response = -1; + } + } + + // Clear the vector of keypoints + kpts.clear(); + + for (size_t i = 0; i < kpts_.size(); i++) { + if (kpts_[i].response != -1) { + kpts.push_back(kpts_[i]); + } + } + + t2 = getTickCount(); + tsubpixel_ = 1000.0*(t2-t1) / getTickFrequency(); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method performs feature suppression based on 2D distance + * @param kpts Vector of keypoints + * @param mdist Maximum distance in pixels +*/ +void KAZE::Feature_Suppression_Distance(std::vector& kpts, const float& mdist) { + + vector aux; + vector to_delete; + float dist = 0.0, x1 = 0.0, y1 = 0.0, x2 = 0.0, y2 = 0.0; + bool found = false; + + for (size_t i = 0; i < kpts.size(); i++) { + x1 = kpts[i].pt.x; + y1 = kpts[i].pt.y; + + for (size_t j = i+1; j < kpts.size(); j++) { + x2 = kpts[j].pt.x; + y2 = kpts[j].pt.y; + dist = sqrt(pow(x1-x2,2)+pow(y1-y2,2)); + + if (dist < mdist) { + if (fabs(kpts[i].response) >= fabs(kpts[j].response)) { + to_delete.push_back(j); + } + else { + to_delete.push_back(i); + break; + } + } + } + } + + for (size_t i = 0; i < kpts.size(); i++) { + found = false; + + for (size_t j = 0; j < to_delete.size(); j++) { + if(i == (size_t)(to_delete[j])) { + found = true; + break; + } + } + + if (found == false) { + aux.push_back(kpts[i]); + } + } + + kpts.clear(); + kpts = aux; + aux.clear(); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the set of descriptors through the nonlinear scale space + * @param kpts Vector of keypoints + * @param desc Matrix with the feature descriptors +*/ +void KAZE::Feature_Description(std::vector &kpts, cv::Mat &desc) { + + double t2 = 0.0, t1 = 0.0; + t1 = getTickCount(); + + // Allocate memory for the matrix of descriptors + if (use_extended_ == true) { + desc = Mat::zeros(kpts.size(),128,CV_32FC1); + } + else { + desc = Mat::zeros(kpts.size(),64,CV_32FC1); + } + + if (use_upright_ == true) { + if (use_extended_ == false) { + if (descriptor_mode_ == 0) { +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (size_t i = 0; i < kpts.size(); i++) { + kpts[i].angle = 0.0; + Get_SURF_Upright_Descriptor_64(kpts[i],desc.ptr(i)); + } + } + else if (descriptor_mode_ == 1) { +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (size_t i = 0; i < kpts.size(); i++) { + kpts[i].angle = 0.0; + Get_MSURF_Upright_Descriptor_64(kpts[i],desc.ptr(i)); + } + } + else if (descriptor_mode_ == 2) { +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (size_t i = 0; i < kpts.size(); i++) { + kpts[i].angle = 0.0; + Get_GSURF_Upright_Descriptor_64(kpts[i],desc.ptr(i)); + } + } + } + else + { + if (descriptor_mode_ == 0) { +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (size_t i = 0; i < kpts.size(); i++ ) { + kpts[i].angle = 0.0; + Get_SURF_Upright_Descriptor_128(kpts[i],desc.ptr(i)); + } + } + else if (descriptor_mode_ == 1) { +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (size_t i = 0; i < kpts.size(); i++) { + kpts[i].angle = 0.0; + Get_MSURF_Upright_Descriptor_128(kpts[i],desc.ptr(i)); + } + } + else if (descriptor_mode_ == 2) { +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (size_t i = 0; i < kpts.size(); i++) { + kpts[i].angle = 0.0; + Get_GSURF_Upright_Descriptor_128(kpts[i],desc.ptr(i)); + } + } + } + } + else { + if (use_extended_ == false) { + if (descriptor_mode_ == 0) { +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (size_t i = 0; i < kpts.size(); i++) { + Compute_Main_Orientation_SURF(kpts[i]); + Get_SURF_Descriptor_64(kpts[i],desc.ptr(i)); + } + } + else if (descriptor_mode_ == 1) { +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (size_t i = 0; i < kpts.size(); i++) { + Compute_Main_Orientation_SURF(kpts[i]); + Get_MSURF_Descriptor_64(kpts[i],desc.ptr(i)); + } + } + else if (descriptor_mode_ == 2) { +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (size_t i = 0; i < kpts.size(); i++) { + Compute_Main_Orientation_SURF(kpts[i]); + Get_GSURF_Descriptor_64(kpts[i],desc.ptr(i)); + } + } + } + else { + if (descriptor_mode_ == 0) { +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (size_t i = 0; i < kpts.size(); i++) { + Compute_Main_Orientation_SURF(kpts[i]); + Get_SURF_Descriptor_128(kpts[i],desc.ptr(i)); + } + } + else if (descriptor_mode_ == 1) { +#ifdef _OPENMP +#pragma omp parallel for +#endif + for(size_t i = 0; i < kpts.size(); i++) { + Compute_Main_Orientation_SURF(kpts[i]); + Get_MSURF_Descriptor_128(kpts[i],desc.ptr(i)); + } + } + else if (descriptor_mode_ == 2) { +#ifdef _OPENMP +#pragma omp parallel for +#endif + for( size_t i = 0; i < kpts.size(); i++ ) { + Compute_Main_Orientation_SURF(kpts[i]); + Get_GSURF_Descriptor_128(kpts[i],desc.ptr(i)); + } + } + } + } + + t2 = getTickCount(); + tdescriptor_ = 1000.0*(t2-t1) / getTickFrequency(); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the main orientation for a given keypoint + * @param kpt Input keypoint + * @note The orientation is computed using a similar approach as described in the + * original SURF method. See Bay et al., Speeded Up Robust Features, ECCV 2006 +*/ +void KAZE::Compute_Main_Orientation_SURF(cv::KeyPoint &kpt) +{ + int ix = 0, iy = 0, idx = 0, s = 0, level = 0; + float xf = 0.0, yf = 0.0, gweight = 0.0; + vector resX(109), resY(109), Ang(109); + + // Variables for computing the dominant direction + float sumX = 0.0, sumY = 0.0, max = 0.0, ang1 = 0.0, ang2 = 0.0; + + // Get the information from the keypoint + xf = kpt.pt.x; + yf = kpt.pt.y; + level = kpt.class_id; + s = fRound(kpt.size/2.0); + + // Calculate derivatives responses for points within radius of 6*scale + for (int i = -6; i <= 6; ++i) { + for (int j = -6; j <= 6; ++j) { + if (i*i + j*j < 36) { + iy = fRound(yf + j*s); + ix = fRound(xf + i*s); + + if (iy >= 0 && iy < img_height_ && ix >= 0 && ix < img_width_ ) { + gweight = gaussian(iy-yf,ix-xf,2.5*s); + resX[idx] = gweight*(*(evolution_[level].Lx.ptr(iy)+ix)); + resY[idx] = gweight*(*(evolution_[level].Ly.ptr(iy)+ix)); + } + else { + resX[idx] = 0.0; + resY[idx] = 0.0; + } + + Ang[idx] = getAngle(resX[idx],resY[idx]); + ++idx; + } + } + } + + // Loop slides pi/3 window around feature point + for (ang1 = 0; ang1 < 2.0*CV_PI; ang1+=0.15f) { + ang2 =(ang1+CV_PI/3.0f > 2.0*CV_PI ? ang1-5.0f*CV_PI/3.0f : ang1+CV_PI/3.0f); + sumX = sumY = 0.f; + + for (size_t k = 0; k < Ang.size(); ++k) { + // Get angle from the x-axis of the sample point + const float & ang = Ang[k]; + + // Determine whether the point is within the window + if (ang1 < ang2 && ang1 < ang && ang < ang2) { + sumX+=resX[k]; + sumY+=resY[k]; + } + else if (ang2 < ang1 && + ((ang > 0 && ang < ang2) || (ang > ang1 && ang < 2.0*CV_PI))) { + sumX+=resX[k]; + sumY+=resY[k]; + } + } + + // if the vector produced from this window is longer than all + // previous vectors then this forms the new dominant direction + if (sumX*sumX + sumY*sumY > max) { + // store largest orientation + max = sumX*sumX + sumY*sumY; + kpt.angle = getAngle(sumX, sumY); + } + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the upright descriptor (no rotation invariant) + * of the provided keypoint + * @param kpt Input keypoint + * @param desc Descriptor vector + * @note Rectangular grid of 20 s x 20 s. Descriptor Length 64. No additional + * Gaussian weighting is performed. The descriptor is inspired from Bay et al., + * Speeded Up Robust Features, ECCV, 2006 +*/ +void KAZE::Get_SURF_Upright_Descriptor_64(const cv::KeyPoint &kpt, float *desc) +{ + float dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0; + float rx = 0.0, ry = 0.0, len = 0.0, xf = 0.0, yf = 0.0, sample_x = 0.0, sample_y = 0.0; + float fx = 0.0, fy = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0, res4 = 0.0; + int x1 = 0, y1 = 0, x2 = 0, y2 = 0, sample_step = 0, pattern_size = 0, dcount = 0; + int dsize = 0, scale = 0, level = 0; + + // Set the descriptor size and the sample and pattern sizes + dsize = 64; + sample_step = 5; + pattern_size = 10; + + // Get the information from the keypoint + yf = kpt.pt.y; + xf = kpt.pt.x; + level = kpt.class_id; + scale = fRound(kpt.size/2.0); + + // Calculate descriptor for this interest point + for (int i = -pattern_size; i < pattern_size; i+=sample_step) { + for (int j = -pattern_size; j < pattern_size; j+=sample_step) { + + dx=dy=mdx=mdy=0.0; + + for (int k = i; k < i + sample_step; k++) { + for (int l = j; l < j + sample_step; l++) { + + sample_y = k*scale + yf; + sample_x = l*scale + xf; + + y1 = (int)(sample_y-.5); + x1 = (int)(sample_x-.5); + + checkDescriptorLimits(x1,y1,img_width_,img_height_); + + y2 = (int)(sample_y+.5); + x2 = (int)(sample_x+.5); + + checkDescriptorLimits(x2,y2,img_width_,img_height_); + + fx = sample_x-x1; + fy = sample_y-y1; + + res1 = *(evolution_[level].Lx.ptr(y1)+x1); + res2 = *(evolution_[level].Lx.ptr(y1)+x2); + res3 = *(evolution_[level].Lx.ptr(y2)+x1); + res4 = *(evolution_[level].Lx.ptr(y2)+x2); + rx = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Ly.ptr(y1)+x1); + res2 = *(evolution_[level].Ly.ptr(y1)+x2); + res3 = *(evolution_[level].Ly.ptr(y2)+x1); + res4 = *(evolution_[level].Ly.ptr(y2)+x2); + ry = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + // Sum the derivatives to the cumulative descriptor + dx += rx; + dy += ry; + mdx += fabs(rx); + mdy += fabs(ry); + } + } + + // Add the values to the descriptor vector + desc[dcount++] = dx; + desc[dcount++] = dy; + desc[dcount++] = mdx; + desc[dcount++] = mdy; + + // Store the current length^2 of the vector + len += dx*dx + dy*dy + mdx*mdx + mdy*mdy; + } + } + + // convert to unit vector + len = sqrt(len); + + for (int i = 0; i < dsize; i++) { + desc[i] /= len; + } + + if (USE_CLIPPING_NORMALIZATION == true) { + clippingDescriptor(desc,dsize,CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO); + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the descriptor of the provided keypoint given the + * main orientation + * @param kpt Input keypoint + * @param desc Descriptor vector + * @note Rectangular grid of 20 s x 20 s. Descriptor Length 64. No additional + * Gaussian weighting is performed. The descriptor is inspired from Bay et al., + * Speeded Up Robust Features, ECCV, 2006 +*/ +void KAZE::Get_SURF_Descriptor_64(const cv::KeyPoint &kpt, float *desc) { + + float dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0; + float rx = 0.0, ry = 0.0, rrx = 0.0, rry = 0.0, len = 0.0, xf = 0.0, yf = 0.0; + float sample_x = 0.0, sample_y = 0.0, co = 0.0, si = 0.0, angle = 0.0; + float fx = 0.0, fy = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0, res4 = 0.0; + int x1 = 0, y1 = 0, x2 = 0, y2 = 0, sample_step = 0, pattern_size = 0, dcount = 0; + int dsize = 0, scale = 0, level = 0; + + // Set the descriptor size and the sample and pattern sizes + dsize = 64; + sample_step = 5; + pattern_size = 10; + + // Get the information from the keypoint + yf = kpt.pt.y; + xf = kpt.pt.x; + scale = fRound(kpt.size/2.0); + angle = kpt.angle; + level = kpt.class_id; + co = cos(angle); + si = sin(angle); + + // Calculate descriptor for this interest point + for (int i = -pattern_size; i < pattern_size; i+=sample_step) { + for (int j = -pattern_size; j < pattern_size; j+=sample_step) { + dx=dy=mdx=mdy=0.0; + + for (int k = i; k < i + sample_step; k++) { + for (int l = j; l < j + sample_step; l++) { + + // Get the coordinates of the sample point on the rotated axis + sample_y = yf + (l*scale*co + k*scale*si); + sample_x = xf + (-l*scale*si + k*scale*co); + + y1 = (int)(sample_y-.5); + x1 = (int)(sample_x-.5); + + checkDescriptorLimits(x1,y1,img_width_,img_height_); + + y2 = (int)(sample_y+.5); + x2 = (int)(sample_x+.5); + + checkDescriptorLimits(x2,y2,img_width_,img_height_); + + fx = sample_x-x1; + fy = sample_y-y1; + + res1 = *(evolution_[level].Lx.ptr(y1)+x1); + res2 = *(evolution_[level].Lx.ptr(y1)+x2); + res3 = *(evolution_[level].Lx.ptr(y2)+x1); + res4 = *(evolution_[level].Lx.ptr(y2)+x2); + rx = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Ly.ptr(y1)+x1); + res2 = *(evolution_[level].Ly.ptr(y1)+x2); + res3 = *(evolution_[level].Ly.ptr(y2)+x1); + res4 = *(evolution_[level].Ly.ptr(y2)+x2); + ry = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + // Get the x and y derivatives on the rotated axis + rry = rx*co + ry*si; + rrx = -rx*si + ry*co; + + // Sum the derivatives to the cumulative descriptor + dx += rrx; + dy += rry; + mdx += fabs(rrx); + mdy += fabs(rry); + } + } + + // Add the values to the descriptor vector + desc[dcount++] = dx; + desc[dcount++] = dy; + desc[dcount++] = mdx; + desc[dcount++] = mdy; + + // Store the current length^2 of the vector + len += dx*dx + dy*dy + mdx*mdx + mdy*mdy; + } + } + + // convert to unit vector + len = sqrt(len); + + for (int i = 0; i < dsize; i++) { + desc[i] /= len; + } + + if (USE_CLIPPING_NORMALIZATION == true) { + clippingDescriptor(desc,dsize,CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO); + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the upright descriptor (not rotation invariant) of + * the provided keypoint + * @param kpt Input keypoint + * @param desc Descriptor vector + * @note Rectangular grid of 24 s x 24 s. Descriptor Length 64. The descriptor is inspired + * from Agrawal et al., CenSurE: Center Surround Extremas for Realtime Feature Detection and Matching, + * ECCV 2008 +*/ +void KAZE::Get_MSURF_Upright_Descriptor_64(const cv::KeyPoint &kpt, float *desc) +{ + float dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0, gauss_s1 = 0.0, gauss_s2 = 0.0; + float rx = 0.0, ry = 0.0, len = 0.0, xf = 0.0, yf = 0.0, ys = 0.0, xs = 0.0; + float sample_x = 0.0, sample_y = 0.0; + int x1 = 0, y1 = 0, sample_step = 0, pattern_size = 0; + int x2 = 0, y2 = 0, kx = 0, ky = 0, i = 0, j = 0, dcount = 0; + float fx = 0.0, fy = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0, res4 = 0.0; + int dsize = 0, scale = 0, level = 0; + + // Subregion centers for the 4x4 gaussian weighting + float cx = -0.5, cy = 0.5; + + // Set the descriptor size and the sample and pattern sizes + dsize = 64; + sample_step = 5; + pattern_size = 12; + + // Get the information from the keypoint + yf = kpt.pt.y; + xf = kpt.pt.x; + scale = fRound(kpt.size/2.0); + level = kpt.class_id; + + i = -8; + + // Calculate descriptor for this interest point + // Area of size 24 s x 24 s + while (i < pattern_size) { + j = -8; + i = i-4; + + cx += 1.0; + cy = -0.5; + + while (j < pattern_size) { + + dx=dy=mdx=mdy=0.0; + cy += 1.0; + j = j-4; + + ky = i + sample_step; + kx = j + sample_step; + + ys = yf + (ky*scale); + xs = xf + (kx*scale); + + for (int k = i; k < i+9; k++) { + for (int l = j; l < j+9; l++) { + + sample_y = k*scale + yf; + sample_x = l*scale + xf; + + //Get the gaussian weighted x and y responses + gauss_s1 = gaussian(xs-sample_x,ys-sample_y,2.5*scale); + + y1 = (int)(sample_y-.5); + x1 = (int)(sample_x-.5); + + checkDescriptorLimits(x1,y1,img_width_,img_height_); + + y2 = (int)(sample_y+.5); + x2 = (int)(sample_x+.5); + + checkDescriptorLimits(x2,y2,img_width_,img_height_); + + fx = sample_x-x1; + fy = sample_y-y1; + + res1 = *(evolution_[level].Lx.ptr(y1)+x1); + res2 = *(evolution_[level].Lx.ptr(y1)+x2); + res3 = *(evolution_[level].Lx.ptr(y2)+x1); + res4 = *(evolution_[level].Lx.ptr(y2)+x2); + rx = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Ly.ptr(y1)+x1); + res2 = *(evolution_[level].Ly.ptr(y1)+x2); + res3 = *(evolution_[level].Ly.ptr(y2)+x1); + res4 = *(evolution_[level].Ly.ptr(y2)+x2); + ry = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + rx = gauss_s1*rx; + ry = gauss_s1*ry; + + // Sum the derivatives to the cumulative descriptor + dx += rx; + dy += ry; + mdx += fabs(rx); + mdy += fabs(ry); + } + } + + // Add the values to the descriptor vector + gauss_s2 = gaussian(cx-2.0f,cy-2.0f,1.5f); + + desc[dcount++] = dx*gauss_s2; + desc[dcount++] = dy*gauss_s2; + desc[dcount++] = mdx*gauss_s2; + desc[dcount++] = mdy*gauss_s2; + + len += (dx*dx + dy*dy + mdx*mdx + mdy*mdy)*gauss_s2*gauss_s2; + + j += 9; + } + + i += 9; + } + + // convert to unit vector + len = sqrt(len); + + for (int i = 0; i < dsize; i++) { + desc[i] /= len; + } + + if (USE_CLIPPING_NORMALIZATION == true) { + clippingDescriptor(desc,dsize,CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO); + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the descriptor of the provided keypoint given the + * main orientation of the keypoint + * @param kpt Input keypoint + * @param desc Descriptor vector + * @note Rectangular grid of 24 s x 24 s. Descriptor Length 64. The descriptor is inspired + * from Agrawal et al., CenSurE: Center Surround Extremas for Realtime Feature Detection and Matching, + * ECCV 2008 +*/ +void KAZE::Get_MSURF_Descriptor_64(const cv::KeyPoint &kpt, float *desc) +{ + float dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0, gauss_s1 = 0.0, gauss_s2 = 0.0; + 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; + float sample_x = 0.0, sample_y = 0.0, co = 0.0, si = 0.0, angle = 0.0; + float fx = 0.0, fy = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0, res4 = 0.0; + int x1 = 0, y1 = 0, x2 = 0, y2 = 0, sample_step = 0, pattern_size = 0; + int kx = 0, ky = 0, i = 0, j = 0, dcount = 0; + int dsize = 0, scale = 0, level = 0; + + // Subregion centers for the 4x4 gaussian weighting + float cx = -0.5, cy = 0.5; + + // Set the descriptor size and the sample and pattern sizes + dsize = 64; + sample_step = 5; + pattern_size = 12; + + // Get the information from the keypoint + yf = kpt.pt.y; + xf = kpt.pt.x; + scale = fRound(kpt.size/2.0); + angle = kpt.angle; + level = kpt.class_id; + co = cos(angle); + si = sin(angle); + + i = -8; + + // Calculate descriptor for this interest point + // Area of size 24 s x 24 s + while (i < pattern_size) { + + j = -8; + i = i-4; + + cx += 1.0; + cy = -0.5; + + while (j < pattern_size) { + + dx=dy=mdx=mdy=0.0; + cy += 1.0; + j = j - 4; + + ky = i + sample_step; + kx = j + sample_step; + + xs = xf + (-kx*scale*si + ky*scale*co); + ys = yf + (kx*scale*co + ky*scale*si); + + for (int k = i; k < i + 9; ++k) { + for (int l = j; l < j + 9; ++l) { + + // Get coords of sample point on the rotated axis + sample_y = yf + (l*scale*co + k*scale*si); + sample_x = xf + (-l*scale*si + k*scale*co); + + // Get the gaussian weighted x and y responses + gauss_s1 = gaussian(xs-sample_x,ys-sample_y,2.5*scale); + y1 = fRound(sample_y-.5); + x1 = fRound(sample_x-.5); + + checkDescriptorLimits(x1,y1,img_width_,img_height_); + + y2 = (int)(sample_y+.5); + x2 = (int)(sample_x+.5); + + checkDescriptorLimits(x2,y2,img_width_,img_height_); + + fx = sample_x-x1; + fy = sample_y-y1; + + res1 = *(evolution_[level].Lx.ptr(y1)+x1); + res2 = *(evolution_[level].Lx.ptr(y1)+x2); + res3 = *(evolution_[level].Lx.ptr(y2)+x1); + res4 = *(evolution_[level].Lx.ptr(y2)+x2); + rx = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Ly.ptr(y1)+x1); + res2 = *(evolution_[level].Ly.ptr(y1)+x2); + res3 = *(evolution_[level].Ly.ptr(y2)+x1); + res4 = *(evolution_[level].Ly.ptr(y2)+x2); + ry = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + // Get the x and y derivatives on the rotated axis + rry = gauss_s1*(rx*co + ry*si); + rrx = gauss_s1*(-rx*si + ry*co); + + // Sum the derivatives to the cumulative descriptor + dx += rrx; + dy += rry; + mdx += fabs(rrx); + mdy += fabs(rry); + } + } + + // Add the values to the descriptor vector + gauss_s2 = gaussian(cx-2.0f,cy-2.0f,1.5f); + desc[dcount++] = dx*gauss_s2; + desc[dcount++] = dy*gauss_s2; + desc[dcount++] = mdx*gauss_s2; + desc[dcount++] = mdy*gauss_s2; + len += (dx*dx + dy*dy + mdx*mdx + mdy*mdy)*gauss_s2*gauss_s2; + j += 9; + } + i += 9; + } + + // convert to unit vector + len = sqrt(len); + + for (int i = 0; i < dsize; i++) { + desc[i] /= len; + } + + if (USE_CLIPPING_NORMALIZATION == true) { + clippingDescriptor(desc,dsize,CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO); + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the upright G-SURF descriptor of the provided keypoint + * given the main orientation + * @param kpt Input keypoint + * @param desc Descriptor vector + * @note Rectangular grid of 20 s x 20 s. Descriptor Length 64. No additional + * G-SURF descriptor as described in Pablo F. Alcantarilla, Luis M. Bergasa and + * Andrew J. Davison, Gauge-SURF Descriptors, Image and Vision Computing 31(1), 2013 +*/ +void KAZE::Get_GSURF_Upright_Descriptor_64(const cv::KeyPoint &kpt, float *desc) +{ + float dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0; + float rx = 0.0, ry = 0.0, rxx = 0.0, rxy = 0.0, ryy = 0.0, len = 0.0, xf = 0.0, yf = 0.0; + float sample_x = 0.0, sample_y = 0.0; + float fx = 0.0, fy = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0, res4 = 0.0; + float lvv = 0.0, lww = 0.0, modg = 0.0; + int x1 = 0, y1 = 0, x2 = 0, y2 = 0, sample_step = 0, pattern_size = 0, dcount = 0; + int dsize = 0, scale = 0, level = 0; + + // Set the descriptor size and the sample and pattern sizes + dsize = 64; + sample_step = 5; + pattern_size = 10; + + // Get the information from the keypoint + yf = kpt.pt.y; + xf = kpt.pt.x; + scale = fRound(kpt.size/2.0); + level = kpt.class_id; + + // Calculate descriptor for this interest point + for (int i = -pattern_size; i < pattern_size; i+=sample_step) { + for (int j = -pattern_size; j < pattern_size; j+=sample_step) { + + dx=dy=mdx=mdy=0.0; + + for (int k = i; k < i + sample_step; k++) { + for (int l = j; l < j + sample_step; l++) { + + // Get the coordinates of the sample point on the rotated axis + sample_y = yf + l*scale; + sample_x = xf + k*scale; + + y1 = (int)(sample_y-.5); + x1 = (int)(sample_x-.5); + + checkDescriptorLimits(x1,y1,img_width_,img_height_); + + y2 = (int)(sample_y+.5); + x2 = (int)(sample_x+.5); + + checkDescriptorLimits(x2,y2,img_width_,img_height_); + + fx = sample_x-x1; + fy = sample_y-y1; + + res1 = *(evolution_[level].Lx.ptr(y1)+x1); + res2 = *(evolution_[level].Lx.ptr(y1)+x2); + res3 = *(evolution_[level].Lx.ptr(y2)+x1); + res4 = *(evolution_[level].Lx.ptr(y2)+x2); + rx = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Ly.ptr(y1)+x1); + res2 = *(evolution_[level].Ly.ptr(y1)+x2); + res3 = *(evolution_[level].Ly.ptr(y2)+x1); + res4 = *(evolution_[level].Ly.ptr(y2)+x2); + ry = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + modg = pow(rx,2) + pow(ry,2); + + if (modg != 0.0) { + + res1 = *(evolution_[level].Lxx.ptr(y1)+x1); + res2 = *(evolution_[level].Lxx.ptr(y1)+x2); + res3 = *(evolution_[level].Lxx.ptr(y2)+x1); + res4 = *(evolution_[level].Lxx.ptr(y2)+x2); + rxx = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Lxy.ptr(y1)+x1); + res2 = *(evolution_[level].Lxy.ptr(y1)+x2); + res3 = *(evolution_[level].Lxy.ptr(y2)+x1); + res4 = *(evolution_[level].Lxy.ptr(y2)+x2); + rxy = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Lyy.ptr(y1)+x1); + res2 = *(evolution_[level].Lyy.ptr(y1)+x2); + res3 = *(evolution_[level].Lyy.ptr(y2)+x1); + res4 = *(evolution_[level].Lyy.ptr(y2)+x2); + ryy = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2) + lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg); + + // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2) + lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg); + } + else { + lww = 0.0; + lvv = 0.0; + } + + // Sum the derivatives to the cumulative descriptor + dx += lww; + dy += lvv; + mdx += fabs(lww); + mdy += fabs(lvv); + } + } + + // Add the values to the descriptor vector + desc[dcount++] = dx; + desc[dcount++] = dy; + desc[dcount++] = mdx; + desc[dcount++] = mdy; + + // Store the current length^2 of the vector + len += dx*dx + dy*dy + mdx*mdx + mdy*mdy; + } + } + + // convert to unit vector + len = sqrt(len); + + for (int i = 0; i < dsize; i++) { + desc[i] /= len; + } + + if (USE_CLIPPING_NORMALIZATION == true) { + clippingDescriptor(desc,dsize,CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO); + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the G-SURF descriptor of the provided keypoint given the + * main orientation + * @param kpt Input keypoint + * @param desc Descriptor vector + * @note Rectangular grid of 20 s x 20 s. Descriptor Length 64. No additional + * G-SURF descriptor as described in Pablo F. Alcantarilla, Luis M. Bergasa and + * Andrew J. Davison, Gauge-SURF Descriptors, Image and Vision Computing 31(1), 2013 +*/ +void KAZE::Get_GSURF_Descriptor_64(const cv::KeyPoint &kpt, float *desc) +{ + float dx = 0.0, dy = 0.0, mdx = 0.0, mdy = 0.0; + float rx = 0.0, ry = 0.0, rxx = 0.0, rxy = 0.0, ryy = 0.0, len = 0.0, xf = 0.0, yf = 0.0; + float sample_x = 0.0, sample_y = 0.0, co = 0.0, si = 0.0, angle = 0.0; + float fx = 0.0, fy = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0, res4 = 0.0; + float lvv = 0.0, lww = 0.0, modg = 0.0; + int x1 = 0, y1 = 0, x2 = 0, y2 = 0, sample_step = 0, pattern_size = 0, dcount = 0; + int dsize = 0, scale = 0, level = 0; + + // Set the descriptor size and the sample and pattern sizes + dsize = 64; + sample_step = 5; + pattern_size = 10; + + // Get the information from the keypoint + yf = kpt.pt.y; + xf = kpt.pt.x; + scale = fRound(kpt.size/2.0); + angle = kpt.angle; + level = kpt.class_id; + co = cos(angle); + si = sin(angle); + + // Calculate descriptor for this interest point + for (int i = -pattern_size; i < pattern_size; i+=sample_step) { + for (int j = -pattern_size; j < pattern_size; j+=sample_step) { + + dx=dy=mdx=mdy=0.0; + + for (int k = i; k < i + sample_step; k++) { + for (int l = j; l < j + sample_step; l++) { + + // Get the coordinates of the sample point on the rotated axis + sample_y = yf + (l*scale*co + k*scale*si); + sample_x = xf + (-l*scale*si + k*scale*co); + + y1 = (int)(sample_y-.5); + x1 = (int)(sample_x-.5); + + checkDescriptorLimits(x1,y1,img_width_,img_height_); + + y2 = (int)(sample_y+.5); + x2 = (int)(sample_x+.5); + + checkDescriptorLimits(x2,y2,img_width_,img_height_); + + fx = sample_x-x1; + fy = sample_y-y1; + + res1 = *(evolution_[level].Lx.ptr(y1)+x1); + res2 = *(evolution_[level].Lx.ptr(y1)+x2); + res3 = *(evolution_[level].Lx.ptr(y2)+x1); + res4 = *(evolution_[level].Lx.ptr(y2)+x2); + rx = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Ly.ptr(y1)+x1); + res2 = *(evolution_[level].Ly.ptr(y1)+x2); + res3 = *(evolution_[level].Ly.ptr(y2)+x1); + res4 = *(evolution_[level].Ly.ptr(y2)+x2); + ry = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + modg = pow(rx,2) + pow(ry,2); + + if (modg != 0.0) { + + res1 = *(evolution_[level].Lxx.ptr(y1)+x1); + res2 = *(evolution_[level].Lxx.ptr(y1)+x2); + res3 = *(evolution_[level].Lxx.ptr(y2)+x1); + res4 = *(evolution_[level].Lxx.ptr(y2)+x2); + rxx = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Lxy.ptr(y1)+x1); + res2 = *(evolution_[level].Lxy.ptr(y1)+x2); + res3 = *(evolution_[level].Lxy.ptr(y2)+x1); + res4 = *(evolution_[level].Lxy.ptr(y2)+x2); + rxy = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Lyy.ptr(y1)+x1); + res2 = *(evolution_[level].Lyy.ptr(y1)+x2); + res3 = *(evolution_[level].Lyy.ptr(y2)+x1); + res4 = *(evolution_[level].Lyy.ptr(y2)+x2); + ryy = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2) + lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg); + + // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2) + lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg); + } + else { + lww = 0.0; + lvv = 0.0; + } + + // Sum the derivatives to the cumulative descriptor + dx += lww; + dy += lvv; + mdx += fabs(lww); + mdy += fabs(lvv); + } + } + + // Add the values to the descriptor vector + desc[dcount++] = dx; + desc[dcount++] = dy; + desc[dcount++] = mdx; + desc[dcount++] = mdy; + + // Store the current length^2 of the vector + len += dx*dx + dy*dy + mdx*mdx + mdy*mdy; + } + } + + // convert to unit vector + len = sqrt(len); + + for (int i = 0; i < dsize; i++) { + desc[i] /= len; + } + + if (USE_CLIPPING_NORMALIZATION == true) { + clippingDescriptor(desc,dsize,CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO); + } + +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the upright extended descriptor (no rotation invariant) + * of the provided keypoint + * @param kpt Input keypoint + * @param desc Descriptor vector + * @note Rectangular grid of 20 s x 20 s. Descriptor Length 128. No additional + * Gaussian weighting is performed. The descriptor is inspired from Bay et al., + * Speeded Up Robust Features, ECCV, 2006 +*/ +void KAZE::Get_SURF_Upright_Descriptor_128(const cv::KeyPoint &kpt, float *desc) +{ + float rx = 0.0, ry = 0.0, len = 0.0, xf = 0.0, yf = 0.0, sample_x = 0.0, sample_y = 0.0; + float fx = 0.0, fy = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0, res4 = 0.0; + float dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0; + float dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0; + int x1 = 0, y1 = 0, x2 = 0, y2 = 0, sample_step = 0, pattern_size = 0, dcount = 0; + int dsize = 0, scale = 0, level = 0; + + // Set the descriptor size and the sample and pattern sizes + dsize = 128; + sample_step = 5; + pattern_size = 10; + + // Get the information from the keypoint + yf = kpt.pt.y; + xf = kpt.pt.x; + scale = fRound(kpt.size/2.0); + level = kpt.class_id; + + // Calculate descriptor for this interest point + for (int i = -pattern_size; i < pattern_size; i+=sample_step) { + for (int j = -pattern_size; j < pattern_size; j+=sample_step) { + + dxp=dxn=mdxp=mdxn=0.0; + dyp=dyn=mdyp=mdyn=0.0; + + for (int k = i; k < i + sample_step; k++) { + for (int l = j; l < j + sample_step; l++) { + + sample_y = k*scale + yf; + sample_x = l*scale + xf; + + y1 = (int)(sample_y-.5); + x1 = (int)(sample_x-.5); + + checkDescriptorLimits(x1,y1,img_width_,img_height_); + + y2 = (int)(sample_y+.5); + x2 = (int)(sample_x+.5); + + checkDescriptorLimits(x2,y2,img_width_,img_height_); + + fx = sample_x-x1; + fy = sample_y-y1; + + res1 = *(evolution_[level].Lx.ptr(y1)+x1); + res2 = *(evolution_[level].Lx.ptr(y1)+x2); + res3 = *(evolution_[level].Lx.ptr(y2)+x1); + res4 = *(evolution_[level].Lx.ptr(y2)+x2); + rx = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Ly.ptr(y1)+x1); + res2 = *(evolution_[level].Ly.ptr(y1)+x2); + res3 = *(evolution_[level].Ly.ptr(y2)+x1); + res4 = *(evolution_[level].Ly.ptr(y2)+x2); + ry = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + // Sum the derivatives to the cumulative descriptor + if (ry >= 0.0) { + dxp += rx; + mdxp += fabs(rx); + } + else { + dxn += rx; + mdxn += fabs(rx); + } + + if (rx >= 0.0) { + dyp += ry; + mdyp += fabs(ry); + } + else { + dyn += ry; + mdyn += fabs(ry); + } + } + } + + // Add the values to the descriptor vector + desc[dcount++] = dxp; + desc[dcount++] = dxn; + desc[dcount++] = mdxp; + desc[dcount++] = mdxn; + desc[dcount++] = dyp; + desc[dcount++] = dyn; + desc[dcount++] = mdyp; + desc[dcount++] = mdyn; + + // Store the current length^2 of the vector + len += dxp*dxp + dxn*dxn + mdxp*mdxp + mdxn*mdxn + + dyp*dyp + dyn*dyn + mdyp*mdyp + mdyn*mdyn; + } + } + + // convert to unit vector + len = sqrt(len); + + for (int i = 0; i < dsize; i++) { + desc[i] /= len; + } + + if (USE_CLIPPING_NORMALIZATION == true) { + clippingDescriptor(desc,dsize,CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO); + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the extended descriptor of the provided keypoint given the + * main orientation + * @param kpt Input keypoint + * @param desc Descriptor vector + * @note Rectangular grid of 20 s x 20 s. Descriptor Length 128. No additional + * Gaussian weighting is performed. The descriptor is inspired from Bay et al., + * Speeded Up Robust Features, ECCV, 2006 +*/ +void KAZE::Get_SURF_Descriptor_128(const cv::KeyPoint &kpt, float *desc) +{ + float rx = 0.0, ry = 0.0, rrx = 0.0, rry = 0.0, len = 0.0, xf = 0.0, yf = 0.0; + float sample_x = 0.0, sample_y = 0.0, co = 0.0, si = 0.0, angle = 0.0; + float fx = 0.0, fy = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0, res4 = 0.0; + float dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0; + float dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0; + int x1 = 0, y1 = 0, x2 = 0, y2 = 0, sample_step = 0, pattern_size = 0, dcount = 0; + int dsize = 0, scale = 0, level = 0; + + // Set the descriptor size and the sample and pattern sizes + dsize = 128; + sample_step = 5; + pattern_size = 10; + + // Get the information from the keypoint + yf = kpt.pt.y; + xf = kpt.pt.x; + scale = fRound(kpt.size/2.0); + angle = kpt.angle; + level = kpt.class_id; + co = cos(angle); + si = sin(angle); + + // Calculate descriptor for this interest point + for (int i = -pattern_size; i < pattern_size; i+=sample_step) { + for (int j = -pattern_size; j < pattern_size; j+=sample_step) { + + dxp=dxn=mdxp=mdxn=0.0; + dyp=dyn=mdyp=mdyn=0.0; + + for (int k = i; k < i + sample_step; k++) { + for (int l = j; l < j + sample_step; l++) { + + // Get the coordinates of the sample point on the rotated axis + sample_y = yf + (l*scale*co + k*scale*si); + sample_x = xf + (-l*scale*si + k*scale*co); + + y1 = (int)(sample_y-.5); + x1 = (int)(sample_x-.5); + + checkDescriptorLimits(x1,y1,img_width_,img_height_); + + y2 = (int)(sample_y+.5); + x2 = (int)(sample_x+.5); + + checkDescriptorLimits(x2,y2,img_width_,img_height_); + + fx = sample_x-x1; + fy = sample_y-y1; + + res1 = *(evolution_[level].Lx.ptr(y1)+x1); + res2 = *(evolution_[level].Lx.ptr(y1)+x2); + res3 = *(evolution_[level].Lx.ptr(y2)+x1); + res4 = *(evolution_[level].Lx.ptr(y2)+x2); + rx = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Ly.ptr(y1)+x1); + res2 = *(evolution_[level].Ly.ptr(y1)+x2); + res3 = *(evolution_[level].Ly.ptr(y2)+x1); + res4 = *(evolution_[level].Ly.ptr(y2)+x2); + ry = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + // Get the x and y derivatives on the rotated axis + rry = rx*co + ry*si; + rrx = -rx*si + ry*co; + + // Sum the derivatives to the cumulative descriptor + if (rry >= 0.0) { + dxp += rrx; + mdxp += fabs(rrx); + } + else { + dxn += rrx; + mdxn += fabs(rrx); + } + + if (rrx >= 0.0) { + dyp += rry; + mdyp += fabs(rry); + } + else { + dyn += rry; + mdyn += fabs(rry); + } + } + } + + // Add the values to the descriptor vector + desc[dcount++] = dxp; + desc[dcount++] = dxn; + desc[dcount++] = mdxp; + desc[dcount++] = mdxn; + desc[dcount++] = dyp; + desc[dcount++] = dyn; + desc[dcount++] = mdyp; + desc[dcount++] = mdyn; + + // Store the current length^2 of the vector + len += dxp*dxp + dxn*dxn + mdxp*mdxp + mdxn*mdxn + + dyp*dyp + dyn*dyn + mdyp*mdyp + mdyn*mdyn; + } + } + + // convert to unit vector + len = sqrt(len); + + for (int i = 0; i < dsize; i++) { + desc[i] /= len; + } + + if (USE_CLIPPING_NORMALIZATION == true) { + clippingDescriptor(desc,dsize,CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO); + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the extended upright descriptor (not rotation invariant) of + * the provided keypoint + * @param kpt Input keypoint + * @param desc Descriptor vector + * @note Rectangular grid of 24 s x 24 s. Descriptor Length 128. The descriptor is inspired + * from Agrawal et al., CenSurE: Center Surround Extremas for Realtime Feature Detection and Matching, + * ECCV 2008 +*/ +void KAZE::Get_MSURF_Upright_Descriptor_128(const cv::KeyPoint &kpt, float *desc) { + + float gauss_s1 = 0.0, gauss_s2 = 0.0; + float rx = 0.0, ry = 0.0, len = 0.0, xf = 0.0, yf = 0.0, ys = 0.0, xs = 0.0; + float sample_x = 0.0, sample_y = 0.0; + int x1 = 0, y1 = 0, sample_step = 0, pattern_size = 0; + int x2 = 0, y2 = 0, kx = 0, ky = 0, i = 0, j = 0, dcount = 0; + float fx = 0.0, fy = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0, res4 = 0.0; + float dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0; + float dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0; + int dsize = 0, scale = 0, level = 0; + + // Subregion centers for the 4x4 gaussian weighting + float cx = -0.5, cy = 0.5; + + // Set the descriptor size and the sample and pattern sizes + dsize = 128; + sample_step = 5; + pattern_size = 12; + + // Get the information from the keypoint + yf = kpt.pt.y; + xf = kpt.pt.x; + scale = fRound(kpt.size/2.0); + level = kpt.class_id; + + i = -8; + + // Calculate descriptor for this interest point + // Area of size 24 s x 24 s + while (i < pattern_size) { + + j = -8; + i = i-4; + + cx += 1.0; + cy = -0.5; + + while (j < pattern_size) { + + dxp=dxn=mdxp=mdxn=0.0; + dyp=dyn=mdyp=mdyn=0.0; + + cy += 1.0; + j = j-4; + + ky = i + sample_step; + kx = j + sample_step; + + ys = yf + (ky*scale); + xs = xf + (kx*scale); + + for (int k = i; k < i+9; k++) { + for (int l = j; l < j+9; l++) { + + sample_y = k*scale + yf; + sample_x = l*scale + xf; + + //Get the gaussian weighted x and y responses + gauss_s1 = gaussian(xs-sample_x,ys-sample_y,2.50*scale); + + y1 = (int)(sample_y-.5); + x1 = (int)(sample_x-.5); + + checkDescriptorLimits(x1,y1,img_width_,img_height_); + + y2 = (int)(sample_y+.5); + x2 = (int)(sample_x+.5); + + checkDescriptorLimits(x2,y2,img_width_,img_height_); + + fx = sample_x-x1; + fy = sample_y-y1; + + res1 = *(evolution_[level].Lx.ptr(y1)+x1); + res2 = *(evolution_[level].Lx.ptr(y1)+x2); + res3 = *(evolution_[level].Lx.ptr(y2)+x1); + res4 = *(evolution_[level].Lx.ptr(y2)+x2); + rx = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Ly.ptr(y1)+x1); + res2 = *(evolution_[level].Ly.ptr(y1)+x2); + res3 = *(evolution_[level].Ly.ptr(y2)+x1); + res4 = *(evolution_[level].Ly.ptr(y2)+x2); + ry = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + rx = gauss_s1*rx; + ry = gauss_s1*ry; + + // Sum the derivatives to the cumulative descriptor + if (ry >= 0.0) { + dxp += rx; + mdxp += fabs(rx); + } + else { + dxn += rx; + mdxn += fabs(rx); + } + + if (rx >= 0.0) { + dyp += ry; + mdyp += fabs(ry); + } + else { + dyn += ry; + mdyn += fabs(ry); + } + } + } + + // Add the values to the descriptor vector + gauss_s2 = gaussian(cx-2.0f,cy-2.0f,1.5f); + + desc[dcount++] = dxp*gauss_s2; + desc[dcount++] = dxn*gauss_s2; + desc[dcount++] = mdxp*gauss_s2; + desc[dcount++] = mdxn*gauss_s2; + desc[dcount++] = dyp*gauss_s2; + desc[dcount++] = dyn*gauss_s2; + desc[dcount++] = mdyp*gauss_s2; + desc[dcount++] = mdyn*gauss_s2; + + // Store the current length^2 of the vector + len += (dxp*dxp + dxn*dxn + mdxp*mdxp + mdxn*mdxn + + dyp*dyp + dyn*dyn + mdyp*mdyp + mdyn*mdyn)*gauss_s2*gauss_s2; + + j += 9; + } + + i += 9; + } + + // convert to unit vector + len = sqrt(len); + + for (int i = 0; i < dsize; i++) { + desc[i] /= len; + } + + if (USE_CLIPPING_NORMALIZATION == true) { + clippingDescriptor(desc,dsize,CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO); + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the extended G-SURF descriptor of the provided keypoint + * given the main orientation of the keypoint + * @param kpt Input keypoint + * @param desc Descriptor vector + * @note Rectangular grid of 24 s x 24 s. Descriptor Length 128. The descriptor is inspired + * from Agrawal et al., CenSurE: Center Surround Extremas for Realtime Feature Detection and Matching, + * ECCV 2008 +*/ +void KAZE::Get_MSURF_Descriptor_128(const cv::KeyPoint &kpt, float *desc) { + + float gauss_s1 = 0.0, gauss_s2 = 0.0; + 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; + float sample_x = 0.0, sample_y = 0.0, co = 0.0, si = 0.0, angle = 0.0; + float fx = 0.0, fy = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0, res4 = 0.0; + float dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0; + float dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0; + int x1 = 0, y1 = 0, x2 = 0, y2 = 0, sample_step = 0, pattern_size = 0; + int kx = 0, ky = 0, i = 0, j = 0, dcount = 0; + int dsize = 0, scale = 0, level = 0; + + // Subregion centers for the 4x4 gaussian weighting + float cx = -0.5, cy = 0.5; + + // Set the descriptor size and the sample and pattern sizes + dsize = 128; + sample_step = 5; + pattern_size = 12; + + // Get the information from the keypoint + yf = kpt.pt.y; + xf = kpt.pt.x; + scale = fRound(kpt.size/2.0); + angle = kpt.angle; + level = kpt.class_id; + co = cos(angle); + si = sin(angle); + + i = -8; + + // Calculate descriptor for this interest point + // Area of size 24 s x 24 s + while (i < pattern_size) { + + j = -8; + i = i-4; + + cx += 1.0; + cy = -0.5; + + while (j < pattern_size) { + + dxp=dxn=mdxp=mdxn=0.0; + dyp=dyn=mdyp=mdyn=0.0; + + cy += 1.0f; + j = j - 4; + + ky = i + sample_step; + kx = j + sample_step; + + xs = xf + (-kx*scale*si + ky*scale*co); + ys = yf + (kx*scale*co + ky*scale*si); + + for (int k = i; k < i + 9; ++k) { + for (int l = j; l < j + 9; ++l) { + + // Get coords of sample point on the rotated axis + sample_y = yf + (l*scale*co + k*scale*si); + sample_x = xf + (-l*scale*si + k*scale*co); + + // Get the gaussian weighted x and y responses + gauss_s1 = gaussian(xs-sample_x,ys-sample_y,2.5*scale); + + y1 = fRound(sample_y-.5); + x1 = fRound(sample_x-.5); + + checkDescriptorLimits(x1,y1,img_width_,img_height_); + + y2 = (int)(sample_y+.5); + x2 = (int)(sample_x+.5); + + checkDescriptorLimits(x2,y2,img_width_,img_height_); + + fx = sample_x-x1; + fy = sample_y-y1; + + res1 = *(evolution_[level].Lx.ptr(y1)+x1); + res2 = *(evolution_[level].Lx.ptr(y1)+x2); + res3 = *(evolution_[level].Lx.ptr(y2)+x1); + res4 = *(evolution_[level].Lx.ptr(y2)+x2); + rx = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Ly.ptr(y1)+x1); + res2 = *(evolution_[level].Ly.ptr(y1)+x2); + res3 = *(evolution_[level].Ly.ptr(y2)+x1); + res4 = *(evolution_[level].Ly.ptr(y2)+x2); + ry = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + // Get the x and y derivatives on the rotated axis + rry = gauss_s1*(rx*co + ry*si); + rrx = gauss_s1*(-rx*si + ry*co); + + // Sum the derivatives to the cumulative descriptor + if (rry >= 0.0) { + dxp += rrx; + mdxp += fabs(rrx); + } + else { + dxn += rrx; + mdxn += fabs(rrx); + } + + if (rrx >= 0.0) { + dyp += rry; + mdyp += fabs(rry); + } + else { + dyn += rry; + mdyn += fabs(rry); + } + } + } + + // Add the values to the descriptor vector + gauss_s2 = gaussian(cx-2.0f,cy-2.0f,1.5f); + + desc[dcount++] = dxp*gauss_s2; + desc[dcount++] = dxn*gauss_s2; + desc[dcount++] = mdxp*gauss_s2; + desc[dcount++] = mdxn*gauss_s2; + desc[dcount++] = dyp*gauss_s2; + desc[dcount++] = dyn*gauss_s2; + desc[dcount++] = mdyp*gauss_s2; + desc[dcount++] = mdyn*gauss_s2; + + // Store the current length^2 of the vector + len += (dxp*dxp + dxn*dxn + mdxp*mdxp + mdxn*mdxn + + dyp*dyp + dyn*dyn + mdyp*mdyp + mdyn*mdyn)*gauss_s2*gauss_s2; + + j += 9; + } + + i += 9; + } + + // convert to unit vector + len = sqrt(len); + + for (int i = 0; i < dsize; i++) { + desc[i] /= len; + } + + if (USE_CLIPPING_NORMALIZATION == true) { + clippingDescriptor(desc,dsize,CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO); + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the G-SURF upright extended descriptor + * (no rotation invariant) of the provided keypoint + * @param kpt Input keypoint + * @param desc Descriptor vector + * @note Rectangular grid of 20 s x 20 s. Descriptor Length 128. No additional + * G-SURF descriptor as described in Pablo F. Alcantarilla, Luis M. Bergasa and + * Andrew J. Davison, Gauge-SURF Descriptors, Image and Vision Computing 31(1), 2013 +*/ +void KAZE::Get_GSURF_Upright_Descriptor_128(const cv::KeyPoint &kpt, float *desc) +{ + float len = 0.0, xf = 0.0, yf = 0.0, sample_x = 0.0, sample_y = 0.0; + float rx = 0.0, ry = 0.0, rxx = 0.0, rxy = 0.0, ryy = 0.0, modg = 0.0; + float fx = 0.0, fy = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0, res4 = 0.0; + float dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0; + float dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0, lvv = 0.0, lww = 0.0; + int x1 = 0, y1 = 0, x2 = 0, y2 = 0, sample_step = 0, pattern_size = 0, dcount = 0; + int dsize = 0, scale = 0, level = 0; + + // Set the descriptor size and the sample and pattern sizes + dsize = 128; + sample_step = 5; + pattern_size = 10; + + // Get the information from the keypoint + yf = kpt.pt.y; + xf = kpt.pt.x; + scale = fRound(kpt.size/2.0); + level = kpt.class_id; + + // Calculate descriptor for this interest point + for (int i = -pattern_size; i < pattern_size; i+=sample_step) { + for(int j = -pattern_size; j < pattern_size; j+=sample_step) { + + dxp=dxn=mdxp=mdxn=0.0; + dyp=dyn=mdyp=mdyn=0.0; + + for (int k = i; k < i + sample_step; k++) { + for (int l = j; l < j + sample_step; l++) { + + sample_y = k*scale + yf; + sample_x = l*scale + xf; + + y1 = (int)(sample_y-.5); + x1 = (int)(sample_x-.5); + + checkDescriptorLimits(x1,y1,img_width_,img_height_); + + y2 = (int)(sample_y+.5); + x2 = (int)(sample_x+.5); + + checkDescriptorLimits(x2,y2,img_width_,img_height_); + + fx = sample_x-x1; + fy = sample_y-y1; + + res1 = *(evolution_[level].Lx.ptr(y1)+x1); + res2 = *(evolution_[level].Lx.ptr(y1)+x2); + res3 = *(evolution_[level].Lx.ptr(y2)+x1); + res4 = *(evolution_[level].Lx.ptr(y2)+x2); + rx = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Ly.ptr(y1)+x1); + res2 = *(evolution_[level].Ly.ptr(y1)+x2); + res3 = *(evolution_[level].Ly.ptr(y2)+x1); + res4 = *(evolution_[level].Ly.ptr(y2)+x2); + ry = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + modg = pow(rx,2) + pow(ry,2); + + if (modg != 0.0) { + + res1 = *(evolution_[level].Lxx.ptr(y1)+x1); + res2 = *(evolution_[level].Lxx.ptr(y1)+x2); + res3 = *(evolution_[level].Lxx.ptr(y2)+x1); + res4 = *(evolution_[level].Lxx.ptr(y2)+x2); + rxx = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Lxy.ptr(y1)+x1); + res2 = *(evolution_[level].Lxy.ptr(y1)+x2); + res3 = *(evolution_[level].Lxy.ptr(y2)+x1); + res4 = *(evolution_[level].Lxy.ptr(y2)+x2); + rxy = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Lyy.ptr(y1)+x1); + res2 = *(evolution_[level].Lyy.ptr(y1)+x2); + res3 = *(evolution_[level].Lyy.ptr(y2)+x1); + res4 = *(evolution_[level].Lyy.ptr(y2)+x2); + ryy = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2) + lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg); + + // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2) + lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg); + } + else { + lww = 0.0; + lvv = 0.0; + } + + // Sum the derivatives to the cumulative descriptor + if (lww >= 0.0) { + dxp += lvv; + mdxp += fabs(lvv); + } + else { + dxn += lvv; + mdxn += fabs(lvv); + } + + if (lvv >= 0.0) { + dyp += lww; + mdyp += fabs(lww); + } + else { + dyn += lww; + mdyn += fabs(lww); + } + } + } + + // Add the values to the descriptor vector + desc[dcount++] = dxp; + desc[dcount++] = dxn; + desc[dcount++] = mdxp; + desc[dcount++] = mdxn; + desc[dcount++] = dyp; + desc[dcount++] = dyn; + desc[dcount++] = mdyp; + desc[dcount++] = mdyn; + + // Store the current length^2 of the vector + len += dxp*dxp + dxn*dxn + mdxp*mdxp + mdxn*mdxn + + dyp*dyp + dyn*dyn + mdyp*mdyp + mdyn*mdyn; + } + } + + // convert to unit vector + len = sqrt(len); + + for (int i = 0; i < dsize; i++) { + desc[i] /= len; + } + + if (USE_CLIPPING_NORMALIZATION == true) { + clippingDescriptor(desc,dsize,CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO); + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method computes the extended descriptor of the provided keypoint given the + * main orientation + * @param kpt Input keypoint + * @param desc Descriptor vector + * @note Rectangular grid of 20 s x 20 s. Descriptor Length 128. No additional + * G-SURF descriptor as described in Pablo F. Alcantarilla, Luis M. Bergasa and + * Andrew J. Davison, Gauge-SURF Descriptors, Image and Vision Computing 31(1), 2013 +*/ +void KAZE::Get_GSURF_Descriptor_128(const cv::KeyPoint &kpt, float *desc) { + + float len = 0.0, xf = 0.0, yf = 0.0; + float rx = 0.0, ry = 0.0, rxx = 0.0, rxy = 0.0, ryy = 0.0; + float sample_x = 0.0, sample_y = 0.0, co = 0.0, si = 0.0, angle = 0.0; + float fx = 0.0, fy = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0, res4 = 0.0; + float dxp = 0.0, dyp = 0.0, mdxp = 0.0, mdyp = 0.0; + float dxn = 0.0, dyn = 0.0, mdxn = 0.0, mdyn = 0.0; + float lvv = 0.0, lww = 0.0, modg = 0.0; + int x1 = 0, y1 = 0, x2 = 0, y2 = 0, sample_step = 0, pattern_size = 0, dcount = 0; + int dsize = 0, scale = 0, level = 0; + + // Set the descriptor size and the sample and pattern sizes + dsize = 128; + sample_step = 5; + pattern_size = 10; + + // Get the information from the keypoint + yf = kpt.pt.y; + xf = kpt.pt.x; + scale = fRound(kpt.size/2.0); + angle = kpt.angle; + level = kpt.class_id; + co = cos(angle); + si = sin(angle); + + // Calculate descriptor for this interest point + for (int i = -pattern_size; i < pattern_size; i+=sample_step) { + for (int j = -pattern_size; j < pattern_size; j+=sample_step) { + + dxp=dxn=mdxp=mdxn=0.0; + dyp=dyn=mdyp=mdyn=0.0; + + for (int k = i; k < i + sample_step; k++) { + for (int l = j; l < j + sample_step; l++) { + + // Get the coordinates of the sample point on the rotated axis + sample_y = yf + (l*scale*co + k*scale*si); + sample_x = xf + (-l*scale*si + k*scale*co); + + y1 = (int)(sample_y-.5); + x1 = (int)(sample_x-.5); + + checkDescriptorLimits(x1,y1,img_width_,img_height_); + + y2 = (int)(sample_y+.5); + x2 = (int)(sample_x+.5); + + checkDescriptorLimits(x2,y2,img_width_,img_height_); + + fx = sample_x-x1; + fy = sample_y-y1; + + res1 = *(evolution_[level].Lx.ptr(y1)+x1); + res2 = *(evolution_[level].Lx.ptr(y1)+x2); + res3 = *(evolution_[level].Lx.ptr(y2)+x1); + res4 = *(evolution_[level].Lx.ptr(y2)+x2); + rx = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Ly.ptr(y1)+x1); + res2 = *(evolution_[level].Ly.ptr(y1)+x2); + res3 = *(evolution_[level].Ly.ptr(y2)+x1); + res4 = *(evolution_[level].Ly.ptr(y2)+x2); + ry = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + modg = pow(rx,2) + pow(ry,2); + + if (modg != 0.0) { + res1 = *(evolution_[level].Lxx.ptr(y1)+x1); + res2 = *(evolution_[level].Lxx.ptr(y1)+x2); + res3 = *(evolution_[level].Lxx.ptr(y2)+x1); + res4 = *(evolution_[level].Lxx.ptr(y2)+x2); + rxx = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Lxy.ptr(y1)+x1); + res2 = *(evolution_[level].Lxy.ptr(y1)+x2); + res3 = *(evolution_[level].Lxy.ptr(y2)+x1); + res4 = *(evolution_[level].Lxy.ptr(y2)+x2); + rxy = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + res1 = *(evolution_[level].Lyy.ptr(y1)+x1); + res2 = *(evolution_[level].Lyy.ptr(y1)+x2); + res3 = *(evolution_[level].Lyy.ptr(y2)+x1); + res4 = *(evolution_[level].Lyy.ptr(y2)+x2); + ryy = (1.0-fx)*(1.0-fy)*res1 + fx*(1.0-fy)*res2 + (1.0-fx)*fy*res3 + fx*fy*res4; + + // Lww = (Lx^2 * Lxx + 2*Lx*Lxy*Ly + Ly^2*Lyy) / (Lx^2 + Ly^2) + lww = (pow(rx,2)*rxx + 2.0*rx*rxy*ry + pow(ry,2)*ryy) / (modg); + + // Lvv = (-2*Lx*Lxy*Ly + Lxx*Ly^2 + Lx^2*Lyy) / (Lx^2 + Ly^2) + lvv = (-2.0*rx*rxy*ry + rxx*pow(ry,2) + pow(rx,2)*ryy) /(modg); + } + else { + lww = 0.0; + lvv = 0.0; + } + + // Sum the derivatives to the cumulative descriptor + if (lww >= 0.0) { + dxp += lvv; + mdxp += fabs(lvv); + } + else { + dxn += lvv; + mdxn += fabs(lvv); + } + + if (lvv >= 0.0) { + dyp += lww; + mdyp += fabs(lww); + } + else { + dyn += lww; + mdyn += fabs(lww); + } + } + } + + // Add the values to the descriptor vector + desc[dcount++] = dxp; + desc[dcount++] = dxn; + desc[dcount++] = mdxp; + desc[dcount++] = mdxn; + desc[dcount++] = dyp; + desc[dcount++] = dyn; + desc[dcount++] = mdyp; + desc[dcount++] = mdyn; + + // Store the current length^2 of the vector + len += dxp*dxp + dxn*dxn + mdxp*mdxp + mdxn*mdxn + + dyp*dyp + dyn*dyn + mdyp*mdyp + mdyn*mdyn; + } + } + + // convert to unit vector + len = sqrt(len); + + for (int i = 0; i < dsize; i++) { + desc[i] /= len; + } + + if (USE_CLIPPING_NORMALIZATION == true) { + clippingDescriptor(desc,dsize,CLIPPING_NORMALIZATION_NITER,CLIPPING_NORMALIZATION_RATIO); + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method performs a scalar non-linear diffusion step using AOS schemes + * @param Ld Image at a given evolution step + * @param Ldprev Image at a previous evolution step + * @param c Conductivity image + * @param stepsize Stepsize for the nonlinear diffusion evolution + * @note If c is constant, the diffusion will be linear + * If c is a matrix of the same size as Ld, the diffusion will be nonlinear + * The stepsize can be arbitrarilly large +*/ +void KAZE::AOS_Step_Scalar(cv::Mat &Ld, const cv::Mat &Ldprev, const cv::Mat &c, const float& stepsize) { + +#ifdef _OPENMP +#pragma omp sections + { +#pragma omp section + { + AOS_Rows(Ldprev,c,stepsize); + } +#pragma omp section + { + AOS_Columns(Ldprev,c,stepsize); + } + } +#else + AOS_Rows(Ldprev,c,stepsize); + AOS_Columns(Ldprev,c,stepsize); +#endif + + Ld = 0.5*(Lty_+Ltx_.t()); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method performs performs 1D-AOS for the image rows + * @param Ldprev Image at a previous evolution step + * @param c Conductivity image + * @param stepsize Stepsize for the nonlinear diffusion evolution +*/ +void KAZE::AOS_Rows(const cv::Mat &Ldprev, const cv::Mat &c, const float& stepsize) { + + // Operate on rows + for (int i = 0; i < qr_.rows; i++) { + for (int j = 0; j < qr_.cols; j++) { + *(qr_.ptr(i)+j) = *(c.ptr(i)+j) + *(c.ptr(i+1)+j); + } + } + + for (int j = 0; j < py_.cols; j++) { + *(py_.ptr(0)+j) = *(qr_.ptr(0)+j); + } + + for (int j = 0; j < py_.cols; j++) { + *(py_.ptr(py_.rows-1)+j) = *(qr_.ptr(qr_.rows-1)+j); + } + + for (int i = 1; i < py_.rows-1; i++) { + for (int j = 0; j < py_.cols; j++) { + *(py_.ptr(i)+j) = *(qr_.ptr(i-1)+j) + *(qr_.ptr(i)+j); + } + } + + // a = 1 + t.*p; (p is -1*p) + // b = -t.*q; + ay_ = 1.0 + stepsize*py_; // p is -1*p + by_ = -stepsize*qr_; + + // Do Thomas algorithm to solve the linear system of equations + Thomas(ay_,by_,Ldprev,Lty_); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method performs performs 1D-AOS for the image columns + * @param Ldprev Image at a previous evolution step + * @param c Conductivity image + * @param stepsize Stepsize for the nonlinear diffusion evolution +*/ +void KAZE::AOS_Columns(const cv::Mat &Ldprev, const cv::Mat &c, const float& stepsize) { + + // Operate on columns + for (int j = 0; j < qc_.cols; j++) { + for (int i = 0; i < qc_.rows; i++) { + *(qc_.ptr(i)+j) = *(c.ptr(i)+j) + *(c.ptr(i)+j+1); + } + } + + for (int i = 0; i < px_.rows; i++) { + *(px_.ptr(i)) = *(qc_.ptr(i)); + } + + for (int i = 0; i < px_.rows; i++) { + *(px_.ptr(i)+px_.cols-1) = *(qc_.ptr(i)+qc_.cols-1); + } + + for (int j = 1; j < px_.cols-1; j++) { + for (int i = 0; i < px_.rows; i++) { + *(px_.ptr(i)+j) = *(qc_.ptr(i)+j-1) + *(qc_.ptr(i)+j); + } + } + + // a = 1 + t.*p'; + ax_ = 1.0 + stepsize*px_.t(); + + // b = -t.*q'; + bx_ = -stepsize*qc_.t(); + + // But take care since we need to transpose the solution!! + Mat Ldprevt = Ldprev.t(); + + // Do Thomas algorithm to solve the linear system of equations + Thomas(ax_,bx_,Ldprevt,Ltx_); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This method does the Thomas algorithm for solving a tridiagonal linear system + * @note The matrix A must be strictly diagonally dominant for a stable solution +*/ +void KAZE::Thomas(const cv::Mat &a, const cv::Mat &b, const cv::Mat &Ld, cv::Mat &x) { + + // Auxiliary variables + int n = a.rows; + Mat m = cv::Mat::zeros(a.rows,a.cols,CV_32F); + Mat l = cv::Mat::zeros(b.rows,b.cols,CV_32F); + Mat y = cv::Mat::zeros(Ld.rows,Ld.cols,CV_32F); + + /** A*x = d; */ + /** / a1 b1 0 0 0 ... 0 \ / x1 \ = / d1 \ */ + /** | c1 a2 b2 0 0 ... 0 | | x2 | = | d2 | */ + /** | 0 c2 a3 b3 0 ... 0 | | x3 | = | d3 | */ + /** | : : : : 0 ... 0 | | : | = | : | */ + /** | : : : : 0 cn-1 an | | xn | = | dn | */ + + /** 1. LU decomposition + / L = / 1 \ U = / m1 r1 \ + / | l1 1 | | m2 r2 | + / | l2 1 | | m3 r3 | + / | : : : | | : : : | + / \ ln-1 1 / \ mn / */ + + for (int j = 0; j < m.cols; j++) { + *(m.ptr(0)+j) = *(a.ptr(0)+j); + } + + for (int j = 0; j < y.cols; j++) { + *(y.ptr(0)+j) = *(Ld.ptr(0)+j); + } + + // 1. Forward substitution L*y = d for y + for (int k = 1; k < n; k++) { + for (int j=0; j < l.cols; j++) { + *(l.ptr(k-1)+j) = *(b.ptr(k-1)+j) / *(m.ptr(k-1)+j); + } + + for (int j=0; j < m.cols; j++) { + *(m.ptr(k)+j) = *(a.ptr(k)+j) - *(l.ptr(k-1)+j)*(*(b.ptr(k-1)+j)); + } + + for (int j=0; j < y.cols; j++) { + *(y.ptr(k)+j) = *(Ld.ptr(k)+j) - *(l.ptr(k-1)+j)*(*(y.ptr(k-1)+j)); + } + } + + // 2. Backward substitution U*x = y + for (int j=0; j < y.cols; j++) { + *(x.ptr(n-1)+j) = (*(y.ptr(n-1)+j))/(*(m.ptr(n-1)+j)); + } + + for (int i = n-2; i >= 0; i--) { + for(int j = 0; j < x.cols; j++) { + *(x.ptr(i)+j) = (*(y.ptr(i)+j) - (*(b.ptr(i)+j))*(*(x.ptr(i+1)+j)))/(*(m.ptr(i)+j)); + } + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function computes the angle from the vector given by (X Y). From 0 to 2*Pi +*/ +inline float getAngle(const float& x, const float& y) { + + if (x >= 0 && y >= 0) + { + return atan(y/x); + } + + if (x < 0 && y >= 0) { + return CV_PI - atan(-y/x); + } + + if(x < 0 && y < 0) { + return CV_PI + atan(y/x); + } + + if(x >= 0 && y < 0) { + return 2.0*CV_PI - atan(-y/x); + } + + return 0; +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function performs descriptor clipping + * @param desc_ Pointer to the descriptor vector + * @param dsize Size of the descriptor vector + * @param iter Number of iterations + * @param ratio Clipping ratio +*/ +inline void clippingDescriptor(float *desc, const int& dsize, const int& niter, const float& ratio) { + + float cratio = ratio / sqrt(dsize); + float len = 0.0; + + for (int i = 0; i < niter; i++) { + len = 0.0; + for (int j = 0; j < dsize; j++) { + if (desc[j] > cratio) { + desc[j] = cratio; + } + else if (desc[j] < -cratio) { + desc[j] = -cratio; + } + len += desc[j]*desc[j]; + } + + // Normalize again + len = sqrt(len); + + for (int j = 0; j < dsize; j++) { + desc[j] = desc[j] / len; + } + } +} + +//************************************************************************************** +//************************************************************************************** + +/** + * @brief This function computes the value of a 2D Gaussian function + * @param x X Position + * @param y Y Position + * @param sig Standard Deviation +*/ +inline float gaussian(const float& x, const float& y, const float& sig) { + return exp(-(x*x+y*y)/(2.0f*sig*sig)); +} + +//************************************************************************************** +//************************************************************************************** + +/** + * @brief This function checks descriptor limits + * @param x X Position + * @param y Y Position + * @param width Image width + * @param height Image height +*/ +inline void checkDescriptorLimits(int &x, int &y, const int& width, const int& height) { + + if (x < 0) { + x = 0; + } + + if (y < 0) { + y = 0; + } + + if (x > width-1) { + x = width-1; + } + + if (y > height-1) { + y = height-1; + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This funtion rounds float to nearest integer + * @param flt Input float + * @return dst Nearest integer + */ +inline int fRound(const float& flt) +{ + return (int)(flt+0.5f); +} diff --git a/modules/features2d/src/kaze/KAZE.h b/modules/features2d/src/kaze/KAZE.h new file mode 100755 index 0000000000..9d489c0480 --- /dev/null +++ b/modules/features2d/src/kaze/KAZE.h @@ -0,0 +1,294 @@ + +/** + * @file KAZE.h + * @brief Main program for detecting and computing descriptors in a nonlinear + * scale space + * @date Jan 21, 2012 + * @author Pablo F. Alcantarilla + */ + +#ifndef KAZE_H_ +#define KAZE_H_ + +//************************************************************************************* +//************************************************************************************* + +// Includes +#include "config.h" +#include "nldiffusion_functions.h" +#include "fed.h" +#include "utils.h" + +//************************************************************************************* +//************************************************************************************* + +// KAZE Class Declaration +class KAZE { + +private: + + // Parameters of the Nonlinear diffusion class + float soffset_; // Base scale offset + float sderivatives_; // Standard deviation of the Gaussian for the nonlinear diff. derivatives + int omax_; // Maximum octave level + int nsublevels_; // Number of sublevels per octave level + int img_width_; // Width of the original image + int img_height_; // Height of the original image + bool save_scale_space_; // For saving scale space images + bool verbosity_; // Verbosity level + std::vector evolution_; // Vector of nonlinear diffusion evolution + float kcontrast_; // The contrast parameter for the scalar nonlinear diffusion + float dthreshold_; // Feature detector threshold response + int diffusivity_; // Diffusivity type, 0->PM G1, 1->PM G2, 2-> Weickert + int descriptor_mode_; // Descriptor mode + bool use_fed_; // Set to true in case we want to use FED for the nonlinear diffusion filtering. Set false for using AOS + bool use_upright_; // Set to true in case we want to use the upright version of the descriptors + bool use_extended_; // Set to true in case we want to use the extended version of the descriptors + + // Vector of keypoint vectors for finding extrema in multiple threads + std::vector > kpts_par_; + + // FED parameters + int ncycles_; // Number of cycles + bool reordering_; // Flag for reordering time steps + std::vector > tsteps_; // Vector of FED dynamic time steps + std::vector nsteps_; // Vector of number of steps per cycle + + // Computation times variables in ms + double tkcontrast_; // Kcontrast factor computation + double tnlscale_; // Nonlinear Scale space generation + double tdetector_; // Feature detector + double tmderivatives_; // Multiscale derivatives computation + double tdresponse_; // Detector response computation + double tdescriptor_; // Feature descriptor + double tsubpixel_; // Subpixel refinement + + // Some auxiliary variables used in the AOS step + cv::Mat Ltx_, Lty_, px_, py_, ax_, ay_, bx_, by_, qr_, qc_; + +public: + + // Constructor + KAZE(KAZEOptions& options); + + // Destructor + ~KAZE(void); + + // Public methods for KAZE interface + void Allocate_Memory_Evolution(void); + int Create_Nonlinear_Scale_Space(const cv::Mat& img); + void Feature_Detection(std::vector& kpts); + void Feature_Description(std::vector& kpts, cv::Mat& desc); + + // Methods for saving the scale space set of images and detector responses + void Save_Nonlinear_Scale_Space(void); + void Save_Detector_Responses(void); + void Save_Flow_Responses(void); + +private: + + // Feature Detection Methods + void Compute_KContrast(const cv::Mat& img, const float& kper); + void Compute_Multiscale_Derivatives(void); + void Compute_Detector_Response(void); + void Determinant_Hessian_Parallel(std::vector& kpts); + void Find_Extremum_Threading(const int& level); + void Do_Subpixel_Refinement(std::vector& kpts); + void Feature_Suppression_Distance(std::vector& kpts, const float& mdist); + + // AOS Methods + void AOS_Step_Scalar(cv::Mat &Ld, const cv::Mat &Ldprev, const cv::Mat &c, const float& stepsize); + void AOS_Rows(const cv::Mat &Ldprev, const cv::Mat &c, const float& stepsize); + void AOS_Columns(const cv::Mat &Ldprev, const cv::Mat &c, const float& stepsize); + void Thomas(const cv::Mat &a, const cv::Mat &b, const cv::Mat &Ld, cv::Mat &x); + + // Feature Description methods + void Compute_Main_Orientation_SURF(cv::KeyPoint& kpt); + + // Descriptor Mode -> 0 SURF 64 + void Get_SURF_Upright_Descriptor_64(const cv::KeyPoint& kpt, float* desc); + void Get_SURF_Descriptor_64(const cv::KeyPoint& kpt, float* desc); + + // Descriptor Mode -> 0 SURF 128 + void Get_SURF_Upright_Descriptor_128(const cv::KeyPoint& kpt, float* desc); + void Get_SURF_Descriptor_128(const cv::KeyPoint& kpt, float* desc); + + // Descriptor Mode -> 1 M-SURF 64 + void Get_MSURF_Upright_Descriptor_64(const cv::KeyPoint& kpt, float* desc); + void Get_MSURF_Descriptor_64(const cv::KeyPoint& kpt, float* desc); + + // Descriptor Mode -> 1 M-SURF 128 + void Get_MSURF_Upright_Descriptor_128(const cv::KeyPoint& kpt, float* desc); + void Get_MSURF_Descriptor_128(const cv::KeyPoint& kpt, float *desc); + + // Descriptor Mode -> 2 G-SURF 64 + void Get_GSURF_Upright_Descriptor_64(const cv::KeyPoint& kpt, float* desc); + void Get_GSURF_Descriptor_64(const cv::KeyPoint& kpt, float *desc); + + // Descriptor Mode -> 2 G-SURF 128 + void Get_GSURF_Upright_Descriptor_128(const cv::KeyPoint& kpt, float* desc); + void Get_GSURF_Descriptor_128(const cv::KeyPoint& kpt, float* desc); + +public: + + // Setters + void Set_Scale_Offset(float soffset) { + soffset_ = soffset; + } + + void Set_SDerivatives(float sderivatives) { + sderivatives_ = sderivatives; + } + + void Set_Octave_Max(int omax) { + omax_ = omax; + } + + void Set_NSublevels(int nsublevels) { + nsublevels_ = nsublevels; + } + + void Set_Save_Scale_Space_Flag(bool save_scale_space) { + save_scale_space_ = save_scale_space; + } + + void Set_Image_Width(int img_width) { + img_width_ = img_width; + } + + void Set_Image_Height(int img_height) { + img_height_ = img_height; + } + + void Set_Verbosity_Level(bool verbosity) { + verbosity_ = verbosity; + } + + void Set_KContrast(float kcontrast) { + kcontrast_ = kcontrast; + } + + void Set_Detector_Threshold(float dthreshold) { + dthreshold_ = dthreshold; + } + + void Set_Diffusivity_Type(int diffusivity) { + diffusivity_ = diffusivity; + } + + void Set_Descriptor_Mode(int descriptor_mode) { + descriptor_mode_ = descriptor_mode; + } + + void Set_Use_FED(bool use_fed) { + use_fed_ = use_fed; + } + + void Set_Upright(bool use_upright) { + use_upright_ = use_upright; + } + + void Set_Extended(bool use_extended) { + use_extended_ = use_extended; + } + + // Getters + float Get_Scale_Offset(void) { + return soffset_; + } + + float Get_SDerivatives(void) { + return sderivatives_; + } + + int Get_Octave_Max(void) { + return omax_; + } + + int Get_NSublevels(void) { + return nsublevels_; + } + + bool Get_Save_Scale_Space_Flag(void) { + return save_scale_space_; + } + + int Get_Image_Width(void) { + return img_width_; + } + + int Get_Image_Height(void) { + return img_height_; + } + + bool Get_Verbosity_Level(void) { + return verbosity_; + } + + float Get_KContrast(void) { + return kcontrast_; + } + + float Get_Detector_Threshold(void) { + return dthreshold_; + } + + int Get_Diffusivity_Type(void) { + return diffusivity_; + } + + int Get_Descriptor_Mode(void) { + return descriptor_mode_; + } + + bool Get_Upright(void) { + return use_upright_; + } + + bool Get_Extended(void) { + return use_extended_; + } + + float Get_Time_KContrast(void) { + return tkcontrast_; + } + + float Get_Time_NLScale(void) { + return tnlscale_; + } + + float Get_Time_Detector(void) { + return tdetector_; + } + + float Get_Time_Multiscale_Derivatives(void) { + return tmderivatives_; + } + + float Get_Time_Detector_Response(void) { + return tdresponse_; + } + + float Get_Time_Descriptor(void) { + return tdescriptor_; + } + + float Get_Time_Subpixel(void) { + return tsubpixel_; + } +}; + +//************************************************************************************* +//************************************************************************************* + +// Inline functions +float getAngle(const float& x, const float& y); +float gaussian(const float& x, const float& y, const float& sig); +void checkDescriptorLimits(int &x, int &y, const int& width, const int& height); +void clippingDescriptor(float *desc, const int& dsize, const int& niter, const float& ratio); +int fRound(const float& flt); + +//************************************************************************************* +//************************************************************************************* + +#endif // KAZE_H_ diff --git a/modules/features2d/src/kaze/config.h b/modules/features2d/src/kaze/config.h new file mode 100644 index 0000000000..ffb41ce826 --- /dev/null +++ b/modules/features2d/src/kaze/config.h @@ -0,0 +1,129 @@ + +/** + * @file config.h + * @brief Configuration file + * @date Dec 27, 2011 + * @author Pablo F. Alcantarilla + */ + +#ifndef _CONFIG_H_ +#define _CONFIG_H_ + +//****************************************************************************** +//****************************************************************************** + +// System Includes +#include +#include +#include +#include +#include +#include +#include + +// OpenCV Includes +#include "precomp.hpp" + +// OpenMP Includes +#ifdef _OPENMP +#include +#else +#define omp_get_thread_num() 0 +#endif + +//************************************************************************************* +//************************************************************************************* + +// Some defines +#define NMAX_CHAR 400 + +// Some default options +const float DEFAULT_SCALE_OFFSET = 1.60; // Base scale offset (sigma units) +const float DEFAULT_OCTAVE_MAX = 4.0; // Maximum octave evolution of the image 2^sigma (coarsest scale sigma units) +const int DEFAULT_NSUBLEVELS = 4; // Default number of sublevels per scale level +const float DEFAULT_DETECTOR_THRESHOLD = 0.001; // Detector response threshold to accept point +const float DEFAULT_MIN_DETECTOR_THRESHOLD = 0.00001; // Minimum Detector response threshold to accept point +const int DEFAULT_DESCRIPTOR_MODE = 1; // Descriptor Mode 0->SURF, 1->M-SURF +const bool DEFAULT_USE_FED = true; // 0->AOS, 1->FED +const bool DEFAULT_UPRIGHT = false; // Upright descriptors, not invariant to rotation +const bool DEFAULT_EXTENDED = false; // Extended descriptor, dimension 128 +const bool DEFAULT_SAVE_SCALE_SPACE = false; // For saving the scale space images +const bool DEFAULT_VERBOSITY = false; // Verbosity level (0->no verbosity) +const bool DEFAULT_SHOW_RESULTS = true; // For showing the output image with the detected features plus some ratios +const bool DEFAULT_SAVE_KEYPOINTS = false; // For saving the list of keypoints + +// Some important configuration variables +const float DEFAULT_SIGMA_SMOOTHING_DERIVATIVES = 1.0; +const float DEFAULT_KCONTRAST = .01; +const float KCONTRAST_PERCENTILE = 0.7; +const int KCONTRAST_NBINS = 300; +const bool COMPUTE_KCONTRAST = true; +const int DEFAULT_DIFFUSIVITY_TYPE = 1; // 0 -> PM G1, 1 -> PM G2, 2 -> Weickert +const bool USE_CLIPPING_NORMALIZATION = false; +const float CLIPPING_NORMALIZATION_RATIO = 1.6; +const int CLIPPING_NORMALIZATION_NITER = 5; + +//************************************************************************************* +//************************************************************************************* + +struct KAZEOptions { + + KAZEOptions() { + // Load the default options + soffset = DEFAULT_SCALE_OFFSET; + omax = DEFAULT_OCTAVE_MAX; + nsublevels = DEFAULT_NSUBLEVELS; + dthreshold = DEFAULT_DETECTOR_THRESHOLD; + use_fed = DEFAULT_USE_FED; + upright = DEFAULT_UPRIGHT; + extended = DEFAULT_EXTENDED; + descriptor = DEFAULT_DESCRIPTOR_MODE; + diffusivity = DEFAULT_DIFFUSIVITY_TYPE; + sderivatives = DEFAULT_SIGMA_SMOOTHING_DERIVATIVES; + save_scale_space = DEFAULT_SAVE_SCALE_SPACE; + save_keypoints = DEFAULT_SAVE_KEYPOINTS; + verbosity = DEFAULT_VERBOSITY; + show_results = DEFAULT_SHOW_RESULTS; + } + + float soffset; + int omax; + int nsublevels; + int img_width; + int img_height; + int diffusivity; + float sderivatives; + float dthreshold; + bool use_fed; + bool upright; + bool extended; + int descriptor; + bool save_scale_space; + bool save_keypoints; + bool verbosity; + bool show_results; +}; + +struct TEvolution { + cv::Mat Lx, Ly; // First order spatial derivatives + cv::Mat Lxx, Lxy, Lyy; // Second order spatial derivatives + cv::Mat Lflow; // Diffusivity image + cv::Mat Lt; // Evolution image + cv::Mat Lsmooth; // Smoothed image + cv::Mat Lstep; // Evolution step update + cv::Mat Ldet; // Detector response + float etime; // Evolution time + float esigma; // Evolution sigma. For linear diffusion t = sigma^2 / 2 + float octave; // Image octave + float sublevel; // Image sublevel in each octave + int sigma_size; // Integer esigma. For computing the feature detector responses +}; + +//************************************************************************************* +//************************************************************************************* + +#endif + + + + diff --git a/modules/features2d/src/kaze/fed.cpp b/modules/features2d/src/kaze/fed.cpp new file mode 100644 index 0000000000..0bd228673a --- /dev/null +++ b/modules/features2d/src/kaze/fed.cpp @@ -0,0 +1,192 @@ +//============================================================================= +// +// fed.cpp +// Authors: Pablo F. Alcantarilla (1), Jesus Nuevo (2) +// Institutions: Georgia Institute of Technology (1) +// TrueVision Solutions (2) +// Date: 15/09/2013 +// Email: pablofdezalc@gmail.com +// +// AKAZE Features Copyright 2013, Pablo F. Alcantarilla, Jesus Nuevo +// All Rights Reserved +// See LICENSE for the license information +//============================================================================= + +/** + * @file fed.cpp + * @brief Functions for performing Fast Explicit Diffusion and building the + * nonlinear scale space + * @date Sep 15, 2013 + * @author Pablo F. Alcantarilla, Jesus Nuevo + * @note This code is derived from FED/FJ library from Grewenig et al., + * The FED/FJ library allows solving more advanced problems + * Please look at the following papers for more information about FED: + * [1] S. Grewenig, J. Weickert, C. Schroers, A. Bruhn. Cyclic Schemes for + * PDE-Based Image Analysis. Technical Report No. 327, Department of Mathematics, + * Saarland University, Saarbrücken, Germany, March 2013 + * [2] S. Grewenig, J. Weickert, A. Bruhn. From box filtering to fast explicit diffusion. + * DAGM, 2010 + * +*/ + +#include "fed.h" + +using namespace std; + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function allocates an array of the least number of time steps such + * that a certain stopping time for the whole process can be obtained and fills + * it with the respective FED time step sizes for one cycle + * The function returns the number of time steps per cycle or 0 on failure + * @param T Desired process stopping time + * @param M Desired number of cycles + * @param tau_max Stability limit for the explicit scheme + * @param reordering Reordering flag + * @param tau The vector with the dynamic step sizes + */ +int fed_tau_by_process_time(const float& T, const int& M, const float& tau_max, + const bool& reordering, std::vector& tau) { + // All cycles have the same fraction of the stopping time + return fed_tau_by_cycle_time(T/(float)M,tau_max,reordering,tau); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function allocates an array of the least number of time steps such + * that a certain stopping time for the whole process can be obtained and fills it + * it with the respective FED time step sizes for one cycle + * The function returns the number of time steps per cycle or 0 on failure + * @param t Desired cycle stopping time + * @param tau_max Stability limit for the explicit scheme + * @param reordering Reordering flag + * @param tau The vector with the dynamic step sizes + */ +int fed_tau_by_cycle_time(const float& t, const float& tau_max, + const bool& reordering, std::vector &tau) { + int n = 0; // Number of time steps + float scale = 0.0; // Ratio of t we search to maximal t + + // Compute necessary number of time steps + n = (int)(ceilf(sqrtf(3.0*t/tau_max+0.25f)-0.5f-1.0e-8f)+ 0.5f); + scale = 3.0*t/(tau_max*(float)(n*(n+1))); + + // Call internal FED time step creation routine + return fed_tau_internal(n,scale,tau_max,reordering,tau); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function allocates an array of time steps and fills it with FED + * time step sizes + * The function returns the number of time steps per cycle or 0 on failure + * @param n Number of internal steps + * @param scale Ratio of t we search to maximal t + * @param tau_max Stability limit for the explicit scheme + * @param reordering Reordering flag + * @param tau The vector with the dynamic step sizes + */ +int fed_tau_internal(const int& n, const float& scale, const float& tau_max, + const bool& reordering, std::vector &tau) { + float c = 0.0, d = 0.0; // Time savers + vector tauh; // Helper vector for unsorted taus + + if (n <= 0) { + return 0; + } + + // Allocate memory for the time step size + tau = vector(n); + + if (reordering) { + tauh = vector(n); + } + + // Compute time saver + c = 1.0f / (4.0f * (float)n + 2.0f); + d = scale * tau_max / 2.0f; + + // Set up originally ordered tau vector + for (int k = 0; k < n; ++k) { + float h = cosf(CV_PI * (2.0f * (float)k + 1.0f) * c); + + if (reordering) { + tauh[k] = d / (h * h); + } + else { + tau[k] = d / (h * h); + } + } + + // Permute list of time steps according to chosen reordering function + int kappa = 0, prime = 0; + + if (reordering == true) { + // Choose kappa cycle with k = n/2 + // This is a heuristic. We can use Leja ordering instead!! + kappa = n / 2; + + // Get modulus for permutation + prime = n + 1; + + while (!fed_is_prime_internal(prime)) { + prime++; + } + + // Perform permutation + for (int k = 0, l = 0; l < n; ++k, ++l) { + int index = 0; + while ((index = ((k+1)*kappa) % prime - 1) >= n) { + k++; + } + + tau[l] = tauh[index]; + } + } + + return n; +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function checks if a number is prime or not + * @param number Number to check if it is prime or not + * @return true if the number is prime + */ +bool fed_is_prime_internal(const int& number) { + bool is_prime = false; + + if (number <= 1) { + return false; + } + else if (number == 1 || number == 2 || number == 3 || number == 5 || number == 7) { + return true; + } + else if ((number % 2) == 0 || (number % 3) == 0 || (number % 5) == 0 || (number % 7) == 0) { + return false; + } + else { + is_prime = true; + int upperLimit = sqrt(number+1.0); + int divisor = 11; + + while (divisor <= upperLimit ) { + if (number % divisor == 0) + { + is_prime = false; + } + + divisor +=2; + } + + return is_prime; + } +} diff --git a/modules/features2d/src/kaze/fed.h b/modules/features2d/src/kaze/fed.h new file mode 100644 index 0000000000..d9e8c49924 --- /dev/null +++ b/modules/features2d/src/kaze/fed.h @@ -0,0 +1,30 @@ +#ifndef FED_H +#define FED_H + +//****************************************************************************** +//****************************************************************************** + +// Includes +#include +#include +#include +#include +#include +#include + +//************************************************************************************* +//************************************************************************************* + +// Declaration of functions +int fed_tau_by_process_time(const float& T, const int& M, const float& tau_max, + const bool& reordering, std::vector& tau); +int fed_tau_by_cycle_time(const float& t, const float& tau_max, + const bool& reordering, std::vector &tau) ; +int fed_tau_internal(const int& n, const float& scale, const float& tau_max, + const bool& reordering, std::vector &tau); +bool fed_is_prime_internal(const int& number); + +//************************************************************************************* +//************************************************************************************* + +#endif // FED_H diff --git a/modules/features2d/src/kaze/nldiffusion_functions.cpp b/modules/features2d/src/kaze/nldiffusion_functions.cpp new file mode 100644 index 0000000000..41a7749058 --- /dev/null +++ b/modules/features2d/src/kaze/nldiffusion_functions.cpp @@ -0,0 +1,386 @@ + +//============================================================================= +// +// nldiffusion_functions.cpp +// Author: Pablo F. Alcantarilla +// Institution: University d'Auvergne +// Address: Clermont Ferrand, France +// Date: 27/12/2011 +// Email: pablofdezalc@gmail.com +// +// KAZE Features Copyright 2012, Pablo F. Alcantarilla +// All Rights Reserved +// See LICENSE for the license information +//============================================================================= + +/** + * @file nldiffusion_functions.cpp + * @brief Functions for non-linear diffusion applications: + * 2D Gaussian Derivatives + * Perona and Malik conductivity equations + * Perona and Malik evolution + * @date Dec 27, 2011 + * @author Pablo F. Alcantarilla + */ + +#include "nldiffusion_functions.h" + +// Namespaces +using namespace std; +using namespace cv; + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function smoothes an image with a Gaussian kernel + * @param src Input image + * @param dst Output image + * @param ksize_x Kernel size in X-direction (horizontal) + * @param ksize_y Kernel size in Y-direction (vertical) + * @param sigma Kernel standard deviation + */ +void gaussian_2D_convolution(const cv::Mat& src, cv::Mat& dst, + int ksize_x, int ksize_y, float sigma) { + + size_t ksize_x_ = 0, ksize_y_ = 0; + + // Compute an appropriate kernel size according to the specified sigma + if (sigma > ksize_x || sigma > ksize_y || ksize_x == 0 || ksize_y == 0) { + ksize_x_ = ceil(2.0*(1.0 + (sigma-0.8)/(0.3))); + ksize_y_ = ksize_x_; + } + + // The kernel size must be and odd number + if ((ksize_x_ % 2) == 0) { + ksize_x_ += 1; + } + + if ((ksize_y_ % 2) == 0) { + ksize_y_ += 1; + } + + // Perform the Gaussian Smoothing with border replication + GaussianBlur(src,dst,Size(ksize_x_,ksize_y_),sigma,sigma,cv::BORDER_REPLICATE); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function computes the Perona and Malik conductivity coefficient g1 + * g1 = exp(-|dL|^2/k^2) + * @param Lx First order image derivative in X-direction (horizontal) + * @param Ly First order image derivative in Y-direction (vertical) + * @param dst Output image + * @param k Contrast factor parameter + */ +void pm_g1(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k) { + cv::exp(-(Lx.mul(Lx) + Ly.mul(Ly))/(k*k),dst); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function computes the Perona and Malik conductivity coefficient g2 + * g2 = 1 / (1 + dL^2 / k^2) + * @param Lx First order image derivative in X-direction (horizontal) + * @param Ly First order image derivative in Y-direction (vertical) + * @param dst Output image + * @param k Contrast factor parameter + */ +void pm_g2(const cv::Mat &Lx, const cv::Mat& Ly, cv::Mat& dst, float k) { + dst = 1./(1. + (Lx.mul(Lx) + Ly.mul(Ly))/(k*k)); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function computes Weickert conductivity coefficient g3 + * @param Lx First order image derivative in X-direction (horizontal) + * @param Ly First order image derivative in Y-direction (vertical) + * @param dst Output image + * @param k Contrast factor parameter + * @note For more information check the following paper: J. Weickert + * Applications of nonlinear diffusion in image processing and computer vision, + * Proceedings of Algorithmy 2000 + */ +void weickert_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k) { + Mat modg; + cv::pow((Lx.mul(Lx) + Ly.mul(Ly))/(k*k),4,modg); + cv::exp(-3.315/modg, dst); + dst = 1.0 - dst; +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function computes a good empirical value for the k contrast factor + * given an input image, the percentile (0-1), the gradient scale and the number of + * bins in the histogram + * @param img Input image + * @param perc Percentile of the image gradient histogram (0-1) + * @param gscale Scale for computing the image gradient histogram + * @param nbins Number of histogram bins + * @param ksize_x Kernel size in X-direction (horizontal) for the Gaussian smoothing kernel + * @param ksize_y Kernel size in Y-direction (vertical) for the Gaussian smoothing kernel + * @return k contrast factor + */ +float compute_k_percentile(const cv::Mat& img, float perc, float gscale, + int nbins, int ksize_x, int ksize_y) { + + int nbin = 0, nelements = 0, nthreshold = 0, k = 0; + float kperc = 0.0, modg = 0.0, lx = 0.0, ly = 0.0; + float npoints = 0.0; + float hmax = 0.0; + + // Create the array for the histogram + float *hist = new float[nbins]; + + // Create the matrices + Mat gaussian = Mat::zeros(img.rows,img.cols,CV_32F); + Mat Lx = Mat::zeros(img.rows,img.cols,CV_32F); + Mat Ly = Mat::zeros(img.rows,img.cols,CV_32F); + + // Set the histogram to zero, just in case + for (int i = 0; i < nbins; i++) { + hist[i] = 0.0; + } + + // Perform the Gaussian convolution + gaussian_2D_convolution(img,gaussian,ksize_x,ksize_y,gscale); + + // Compute the Gaussian derivatives Lx and Ly + Scharr(gaussian,Lx,CV_32F,1,0,1,0,cv::BORDER_DEFAULT); + Scharr(gaussian,Ly,CV_32F,0,1,1,0,cv::BORDER_DEFAULT); + + // Skip the borders for computing the histogram + for (int i = 1; i < gaussian.rows-1; i++) { + for (int j = 1; j < gaussian.cols-1; j++) { + lx = *(Lx.ptr(i)+j); + ly = *(Ly.ptr(i)+j); + modg = sqrt(lx*lx + ly*ly); + + // Get the maximum + if (modg > hmax) { + hmax = modg; + } + } + } + + // Skip the borders for computing the histogram + for (int i = 1; i < gaussian.rows-1; i++) { + for (int j = 1; j < gaussian.cols-1; j++) { + lx = *(Lx.ptr(i)+j); + ly = *(Ly.ptr(i)+j); + modg = sqrt(lx*lx + ly*ly); + + // Find the correspondent bin + if (modg != 0.0) { + nbin = floor(nbins*(modg/hmax)); + + if (nbin == nbins) { + nbin--; + } + + hist[nbin]++; + npoints++; + } + } + } + + // Now find the perc of the histogram percentile + nthreshold = (size_t)(npoints*perc); + + + for (k = 0; nelements < nthreshold && k < nbins; k++) { + nelements = nelements + hist[k]; + } + + if (nelements < nthreshold) { + kperc = 0.03; + } + else { + kperc = hmax*((float)(k)/(float)nbins); + } + + delete hist; + return kperc; +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function computes Scharr image derivatives + * @param src Input image + * @param dst Output image + * @param xorder Derivative order in X-direction (horizontal) + * @param yorder Derivative order in Y-direction (vertical) + * @param scale Scale factor or derivative size + */ +void compute_scharr_derivatives(const cv::Mat& src, cv::Mat& dst, + int xorder, int yorder, int scale) { + Mat kx, ky; + compute_derivative_kernels(kx,ky,xorder,yorder,scale); + sepFilter2D(src,dst,CV_32F,kx,ky); +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief Compute derivative kernels for sizes different than 3 + * @param _kx Horizontal kernel values + * @param _ky Vertical kernel values + * @param dx Derivative order in X-direction (horizontal) + * @param dy Derivative order in Y-direction (vertical) + * @param scale_ Scale factor or derivative size + */ +void compute_derivative_kernels(cv::OutputArray _kx, cv::OutputArray _ky, + int dx, int dy, int scale) { + + int ksize = 3 + 2*(scale-1); + + // The standard Scharr kernel + if (scale == 1) { + getDerivKernels(_kx,_ky,dx,dy,0,true,CV_32F); + return; + } + + _kx.create(ksize,1,CV_32F,-1,true); + _ky.create(ksize,1,CV_32F,-1,true); + Mat kx = _kx.getMat(); + Mat ky = _ky.getMat(); + + float w = 10.0/3.0; + float norm = 1.0/(2.0*scale*(w+2.0)); + + for (int k = 0; k < 2; k++) { + Mat* kernel = k == 0 ? &kx : &ky; + int order = k == 0 ? dx : dy; + std::vector kerI(ksize); + + for (int t=0; trows,kernel->cols,CV_32F,&kerI[0]); + temp.copyTo(*kernel); + } +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function performs a scalar non-linear diffusion step + * @param Ld2 Output image in the evolution + * @param c Conductivity image + * @param Lstep Previous image in the evolution + * @param stepsize The step size in time units + * @note Forward Euler Scheme 3x3 stencil + * The function c is a scalar value that depends on the gradient norm + * dL_by_ds = d(c dL_by_dx)_by_dx + d(c dL_by_dy)_by_dy + */ +void nld_step_scalar(cv::Mat& Ld, const cv::Mat& c, cv::Mat& Lstep, float stepsize) { + +#ifdef _OPENMP +#pragma omp parallel for schedule(dynamic) +#endif + for (int i = 1; i < Lstep.rows-1; i++) { + for (int j = 1; j < Lstep.cols-1; j++) { + float xpos = ((*(c.ptr(i)+j))+(*(c.ptr(i)+j+1)))*((*(Ld.ptr(i)+j+1))-(*(Ld.ptr(i)+j))); + float xneg = ((*(c.ptr(i)+j-1))+(*(c.ptr(i)+j)))*((*(Ld.ptr(i)+j))-(*(Ld.ptr(i)+j-1))); + float ypos = ((*(c.ptr(i)+j))+(*(c.ptr(i+1)+j)))*((*(Ld.ptr(i+1)+j))-(*(Ld.ptr(i)+j))); + float yneg = ((*(c.ptr(i-1)+j))+(*(c.ptr(i)+j)))*((*(Ld.ptr(i)+j))-(*(Ld.ptr(i-1)+j))); + *(Lstep.ptr(i)+j) = 0.5*stepsize*(xpos-xneg + ypos-yneg); + } + } + + for (int j = 1; j < Lstep.cols-1; j++) { + float xpos = ((*(c.ptr(0)+j))+(*(c.ptr(0)+j+1)))*((*(Ld.ptr(0)+j+1))-(*(Ld.ptr(0)+j))); + float xneg = ((*(c.ptr(0)+j-1))+(*(c.ptr(0)+j)))*((*(Ld.ptr(0)+j))-(*(Ld.ptr(0)+j-1))); + float ypos = ((*(c.ptr(0)+j))+(*(c.ptr(1)+j)))*((*(Ld.ptr(1)+j))-(*(Ld.ptr(0)+j))); + float yneg = ((*(c.ptr(0)+j))+(*(c.ptr(0)+j)))*((*(Ld.ptr(0)+j))-(*(Ld.ptr(0)+j))); + *(Lstep.ptr(0)+j) = 0.5*stepsize*(xpos-xneg + ypos-yneg); + } + + for (int j = 1; j < Lstep.cols-1; j++) { + float xpos = ((*(c.ptr(Lstep.rows-1)+j))+(*(c.ptr(Lstep.rows-1)+j+1)))*((*(Ld.ptr(Lstep.rows-1)+j+1))-(*(Ld.ptr(Lstep.rows-1)+j))); + float xneg = ((*(c.ptr(Lstep.rows-1)+j-1))+(*(c.ptr(Lstep.rows-1)+j)))*((*(Ld.ptr(Lstep.rows-1)+j))-(*(Ld.ptr(Lstep.rows-1)+j-1))); + float ypos = ((*(c.ptr(Lstep.rows-1)+j))+(*(c.ptr(Lstep.rows-1)+j)))*((*(Ld.ptr(Lstep.rows-1)+j))-(*(Ld.ptr(Lstep.rows-1)+j))); + float yneg = ((*(c.ptr(Lstep.rows-2)+j))+(*(c.ptr(Lstep.rows-1)+j)))*((*(Ld.ptr(Lstep.rows-1)+j))-(*(Ld.ptr(Lstep.rows-2)+j))); + *(Lstep.ptr(Lstep.rows-1)+j) = 0.5*stepsize*(xpos-xneg + ypos-yneg); + } + + for (int i = 1; i < Lstep.rows-1; i++) { + float xpos = ((*(c.ptr(i)))+(*(c.ptr(i)+1)))*((*(Ld.ptr(i)+1))-(*(Ld.ptr(i)))); + float xneg = ((*(c.ptr(i)))+(*(c.ptr(i))))*((*(Ld.ptr(i)))-(*(Ld.ptr(i)))); + float ypos = ((*(c.ptr(i)))+(*(c.ptr(i+1))))*((*(Ld.ptr(i+1)))-(*(Ld.ptr(i)))); + float yneg = ((*(c.ptr(i-1)))+(*(c.ptr(i))))*((*(Ld.ptr(i)))-(*(Ld.ptr(i-1)))); + *(Lstep.ptr(i)) = 0.5*stepsize*(xpos-xneg + ypos-yneg); + } + + for (int i = 1; i < Lstep.rows-1; i++) { + float xpos = ((*(c.ptr(i)+Lstep.cols-1))+(*(c.ptr(i)+Lstep.cols-1)))*((*(Ld.ptr(i)+Lstep.cols-1))-(*(Ld.ptr(i)+Lstep.cols-1))); + float xneg = ((*(c.ptr(i)+Lstep.cols-2))+(*(c.ptr(i)+Lstep.cols-1)))*((*(Ld.ptr(i)+Lstep.cols-1))-(*(Ld.ptr(i)+Lstep.cols-2))); + float ypos = ((*(c.ptr(i)+Lstep.cols-1))+(*(c.ptr(i+1)+Lstep.cols-1)))*((*(Ld.ptr(i+1)+Lstep.cols-1))-(*(Ld.ptr(i)+Lstep.cols-1))); + float yneg = ((*(c.ptr(i-1)+Lstep.cols-1))+(*(c.ptr(i)+Lstep.cols-1)))*((*(Ld.ptr(i)+Lstep.cols-1))-(*(Ld.ptr(i-1)+Lstep.cols-1))); + *(Lstep.ptr(i)+Lstep.cols-1) = 0.5*stepsize*(xpos-xneg + ypos-yneg); + } + + Ld = Ld + Lstep; +} + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function checks if a given pixel is a maximum in a local neighbourhood + * @param img Input image where we will perform the maximum search + * @param dsize Half size of the neighbourhood + * @param value Response value at (x,y) position + * @param row Image row coordinate + * @param col Image column coordinate + * @param same_img Flag to indicate if the image value at (x,y) is in the input image + * @return 1->is maximum, 0->otherwise + */ +bool check_maximum_neighbourhood(const cv::Mat& img, int dsize, float value, + int row, int col, bool same_img) { + + bool response = true; + + for (int i = row-dsize; i <= row+dsize; i++) { + for (int j = col-dsize; j <= col+dsize; j++) { + if (i >= 0 && i < img.rows && j >= 0 && j < img.cols) { + if (same_img == true) { + if (i != row || j != col) { + if ((*(img.ptr(i)+j)) > value) { + response = false; + return response; + } + } + } + else { + if ((*(img.ptr(i)+j)) > value) { + response = false; + return response; + } + } + } + } + } + + return response; +} diff --git a/modules/features2d/src/kaze/nldiffusion_functions.h b/modules/features2d/src/kaze/nldiffusion_functions.h new file mode 100755 index 0000000000..d0ece89571 --- /dev/null +++ b/modules/features2d/src/kaze/nldiffusion_functions.h @@ -0,0 +1,51 @@ + +/** + * @file nldiffusion_functions.h + * @brief Functions for non-linear diffusion applications: + * 2D Gaussian Derivatives + * Perona and Malik conductivity equations + * Perona and Malik evolution + * @date Dec 27, 2011 + * @author Pablo F. Alcantarilla + */ + +#ifndef NLDIFFUSION_FUNCTIONS_H_ +#define NLDIFFUSION_FUNCTIONS_H_ + +//****************************************************************************** +//****************************************************************************** + +// Includes +#include "config.h" + +//************************************************************************************* +//************************************************************************************* + +// Gaussian 2D convolution +void gaussian_2D_convolution(const cv::Mat& src, cv::Mat& dst, + int ksize_x, int ksize_y, float sigma); + +// Diffusivity functions +void pm_g1(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k); +void pm_g2(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k); +void weickert_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k); +float compute_k_percentile(const cv::Mat& img, float perc, float gscale, + int nbins, int ksize_x, int ksize_y); + +// Image derivatives +void compute_scharr_derivatives(const cv::Mat& src, cv::Mat& dst, + int xorder, int yorder, int scale); +void compute_derivative_kernels(cv::OutputArray _kx, cv::OutputArray _ky, + int dx, int dy, int scale); + +// Nonlinear diffusion filtering scalar step +void nld_step_scalar(cv::Mat& Ld, const cv::Mat& c, cv::Mat& Lstep, float stepsize); + +// For non-maxima suppresion +bool check_maximum_neighbourhood(const cv::Mat& img, int dsize, float value, + int row, int col, bool same_img); + +//************************************************************************************* +//************************************************************************************* + +#endif // NLDIFFUSION_FUNCTIONS_H_ diff --git a/modules/features2d/src/kaze/utils.cpp b/modules/features2d/src/kaze/utils.cpp new file mode 100644 index 0000000000..7b55ac45f7 --- /dev/null +++ b/modules/features2d/src/kaze/utils.cpp @@ -0,0 +1,92 @@ + +//============================================================================= +// +// utils.cpp +// Author: Pablo F. Alcantarilla +// Institution: University d'Auvergne +// Address: Clermont Ferrand, France +// Date: 29/12/2011 +// Email: pablofdezalc@gmail.com +// +// KAZE Features Copyright 2012, Pablo F. Alcantarilla +// All Rights Reserved +// See LICENSE for the license information +//============================================================================= + +/** + * @file utils.cpp + * @brief Some useful functions + * @date Dec 29, 2011 + * @author Pablo F. Alcantarilla + */ + +#include "utils.h" + +using namespace std; +using namespace cv; + +//************************************************************************************* +//************************************************************************************* + +/** + * @brief This function copies the input image and converts the scale of the copied + * image prior visualization + * @param src Input image + * @param dst Output image + */ +void copy_and_convert_scale(const cv::Mat& src, cv::Mat& dst) { + + float min_val = 0, max_val = 0; + + src.copyTo(dst); + compute_min_32F(dst,min_val); + + dst = dst - min_val; + + compute_max_32F(dst,max_val); + dst = dst / max_val; +} + +//************************************************************************************* +//************************************************************************************* + +/* +void show_input_options_help(int example) { + + fflush(stdout); + + cout << endl; + cout << endl; + cout << "KAZE Features" << endl; + cout << "***********************************************************" << endl; + cout << "For running the program you need to type in the command line the following arguments: " << endl; + + if (example == 0) { + cout << "./kaze_features img.jpg [options]" << endl; + } + else if (example == 1) { + cout << "./kaze_match img1.jpg img2.pgm homography.txt [options]" << endl; + } + else if (example == 2) { + cout << "./kaze_compare img1.jpg img2.pgm homography.txt [options]" << endl; + } + + cout << endl; + cout << "The options are not mandatory. In case you do not specify additional options, default arguments will be used" << endl << endl; + cout << "Here is a description of the additional options: " << endl; + cout << "--verbose " << "\t\t if verbosity is required" << endl; + cout << "--help" << "\t\t for showing the command line options" << endl; + cout << "--soffset" << "\t\t the base scale offset (sigma units)" << endl; + cout << "--omax" << "\t\t maximum octave evolution of the image 2^sigma (coarsest scale)" << endl; + cout << "--nsublevels" << "\t\t number of sublevels per octave" << endl; + cout << "--dthreshold" << "\t\t Feature detector threshold response for accepting points (0.001 can be a good value)" << endl; + cout << "--descriptor" << "\t\t Descriptor Type 0 -> SURF, 1 -> M-SURF, 2 -> G-SURF" << endl; + cout << "--use_fed" "\t\t 1 -> Use FED, 0 -> Use AOS for the nonlinear diffusion filtering" << endl; + cout << "--upright" << "\t\t 0 -> Rotation Invariant, 1 -> No Rotation Invariant" << endl; + cout << "--extended" << "\t\t 0 -> Normal Descriptor (64), 1 -> Extended Descriptor (128)" << endl; + cout << "--output keypoints.txt" << "\t\t For saving the detected keypoints into a .txt file" << endl; + cout << "--save_scale_space" << "\t\t 1 in case we want to save the nonlinear scale space images. 0 otherwise" << endl; + cout << "--show_results" << "\t\t 1 in case we want to show detection results. 0 otherwise" << endl; + cout << endl; +} +*/ \ No newline at end of file diff --git a/modules/features2d/src/kaze/utils.h b/modules/features2d/src/kaze/utils.h new file mode 100644 index 0000000000..848bfe3f53 --- /dev/null +++ b/modules/features2d/src/kaze/utils.h @@ -0,0 +1,41 @@ + +/** + * @file utils.h + * @brief Some useful functions + * @date Dec 29, 2011 + * @author Pablo F. Alcantarilla + */ + +#ifndef UTILS_H_ +#define UTILS_H_ + +//****************************************************************************** +//****************************************************************************** + +// OPENCV Includes +#include "precomp.hpp" + +// System Includes +#include +#include +#include +#include +#include +#include +#include +#include +#include + +//************************************************************************************* +//************************************************************************************* + +// Declaration of Functions +void compute_min_32F(const cv::Mat& src, float& value); +void compute_max_32F(const cv::Mat& src, float& value); +void convert_scale(cv::Mat& src); +void copy_and_convert_scale(const cv::Mat &src, cv::Mat& dst); + +//************************************************************************************* +//************************************************************************************* + +#endif // UTILS_H_