From 86888bdf7797315731603fd336905793e81120ef Mon Sep 17 00:00:00 2001 From: Ievgen Khvedchenia Date: Thu, 24 Apr 2014 22:01:45 +0100 Subject: [PATCH] Replace swap with clear (more efficient) --- modules/features2d/src/akaze/AKAZE.cpp | 2903 ++++++++++++------------ 1 file changed, 1452 insertions(+), 1451 deletions(-) diff --git a/modules/features2d/src/akaze/AKAZE.cpp b/modules/features2d/src/akaze/AKAZE.cpp index 661a1cad8b..617a16d4ee 100644 --- a/modules/features2d/src/akaze/AKAZE.cpp +++ b/modules/features2d/src/akaze/AKAZE.cpp @@ -18,81 +18,81 @@ using namespace cv; * @brief AKAZEFeatures constructor with input options * @param options AKAZEFeatures configuration options * @note This constructor allocates memory for the nonlinear scale space -*/ + */ AKAZEFeatures::AKAZEFeatures(const AKAZEOptions& options) : options_(options) { - ncycles_ = 0; - reordering_ = true; + ncycles_ = 0; + reordering_ = true; - if (options_.descriptor_size > 0 && options_.descriptor >= MLDB_UPRIGHT) { - generateDescriptorSubsample(descriptorSamples_, descriptorBits_, options_.descriptor_size, - options_.descriptor_pattern_size, options_.descriptor_channels); - } + if (options_.descriptor_size > 0 && options_.descriptor >= MLDB_UPRIGHT) { + generateDescriptorSubsample(descriptorSamples_, descriptorBits_, options_.descriptor_size, + options_.descriptor_pattern_size, options_.descriptor_channels); + } - Allocate_Memory_Evolution(); + Allocate_Memory_Evolution(); } /* ************************************************************************* */ /** * @brief AKAZEFeatures destructor -*/ + */ AKAZEFeatures::~AKAZEFeatures(void) { - evolution_.clear(); + evolution_.clear(); } /* ************************************************************************* */ /** * @brief This method allocates the memory for the nonlinear diffusion evolution -*/ + */ void AKAZEFeatures::Allocate_Memory_Evolution(void) { - float rfactor = 0.0; - int level_height = 0, level_width = 0; + 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 <= options_.omax-1; i++) { - rfactor = 1.0/pow(2.f, i); - level_height = (int)(options_.img_height*rfactor); - level_width = (int)(options_.img_width*rfactor); + // Allocate the dimension of the matrices for the evolution + for (int i = 0; i <= options_.omax - 1; i++) { + rfactor = 1.0 / pow(2.f, i); + level_height = (int)(options_.img_height*rfactor); + level_width = (int)(options_.img_width*rfactor); - // Smallest possible octave and allow one scale if the image is small - if ((level_width < 80 || level_height < 40) && i != 0) { - options_.omax = i; - break; + // Smallest possible octave and allow one scale if the image is small + if ((level_width < 80 || level_height < 40) && i != 0) { + options_.omax = i; + break; + } + + for (int j = 0; j < options_.nsublevels; j++) { + TEvolution step; + step.Lx = cv::Mat::zeros(level_height, level_width, CV_32F); + step.Ly = cv::Mat::zeros(level_height, level_width, CV_32F); + step.Lxx = cv::Mat::zeros(level_height, level_width, CV_32F); + step.Lxy = cv::Mat::zeros(level_height, level_width, CV_32F); + step.Lyy = cv::Mat::zeros(level_height, level_width, CV_32F); + step.Lt = cv::Mat::zeros(level_height, level_width, CV_32F); + step.Ldet = cv::Mat::zeros(level_height, level_width, CV_32F); + step.Lflow = cv::Mat::zeros(level_height, level_width, CV_32F); + step.Lstep = cv::Mat::zeros(level_height, level_width, CV_32F); + step.esigma = options_.soffset*pow(2.f, (float)(j) / (float)(options_.nsublevels) + i); + step.sigma_size = fRound(step.esigma); + step.etime = 0.5*(step.esigma*step.esigma); + step.octave = i; + step.sublevel = j; + evolution_.push_back(step); + } } - for (int j = 0; j < options_.nsublevels; j++) { - TEvolution step; - step.Lx = cv::Mat::zeros(level_height, level_width, CV_32F); - step.Ly = cv::Mat::zeros(level_height, level_width, CV_32F); - step.Lxx = cv::Mat::zeros(level_height, level_width, CV_32F); - step.Lxy = cv::Mat::zeros(level_height, level_width, CV_32F); - step.Lyy = cv::Mat::zeros(level_height, level_width, CV_32F); - step.Lt = cv::Mat::zeros(level_height, level_width, CV_32F); - step.Ldet = cv::Mat::zeros(level_height, level_width, CV_32F); - step.Lflow = cv::Mat::zeros(level_height, level_width, CV_32F); - step.Lstep = cv::Mat::zeros(level_height, level_width, CV_32F); - step.esigma = options_.soffset*pow(2.f, (float)(j)/(float)(options_.nsublevels) + i); - step.sigma_size = fRound(step.esigma); - step.etime = 0.5*(step.esigma*step.esigma); - step.octave = i; - step.sublevel = j; - evolution_.push_back(step); + // Allocate memory for the number of cycles and time steps + 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_++; } - } - - // Allocate memory for the number of cycles and time steps - 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_++; - } } /* ************************************************************************* */ @@ -100,364 +100,365 @@ void AKAZEFeatures::Allocate_Memory_Evolution(void) { * @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 AKAZEFeatures::Create_Nonlinear_Scale_Space(const cv::Mat& img) { - double t1 = 0.0, t2 = 0.0; + double t1 = 0.0, t2 = 0.0; - if (evolution_.size() == 0) { - cerr << "Error generating the nonlinear scale space!!" << endl; - cerr << "Firstly you need to call AKAZEFeatures::Allocate_Memory_Evolution()" << endl; - return -1; - } - - t1 = cv::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, options_.soffset); - evolution_[0].Lt.copyTo(evolution_[0].Lsmooth); - - // First compute the kcontrast factor - options_.kcontrast = compute_k_percentile(img, options_.kcontrast_percentile, - 1.0, options_.kcontrast_nbins, 0, 0); - - t2 = cv::getTickCount(); - timing_.kcontrast = 1000.0*(t2-t1) / cv::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); - options_.kcontrast = options_.kcontrast*0.75; - } - else { - evolution_[i-1].Lt.copyTo(evolution_[i].Lt); + if (evolution_.size() == 0) { + cerr << "Error generating the nonlinear scale space!!" << endl; + cerr << "Firstly you need to call AKAZEFeatures::Allocate_Memory_Evolution()" << endl; + return -1; } - gaussian_2D_convolution(evolution_[i].Lt, evolution_[i].Lsmooth, 0, 0, 1.0); + t1 = cv::getTickCount(); - // 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); + // 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, options_.soffset); + evolution_[0].Lt.copyTo(evolution_[0].Lsmooth); - // Compute the conductivity equation - switch (options_.diffusivity) { - case PM_G1: - pm_g1(evolution_[i].Lx, evolution_[i].Ly, evolution_[i].Lflow, options_.kcontrast); - break; - case PM_G2: - pm_g2(evolution_[i].Lx, evolution_[i].Ly, evolution_[i].Lflow, options_.kcontrast); - break; - case WEICKERT: - weickert_diffusivity(evolution_[i].Lx, evolution_[i].Ly, evolution_[i].Lflow, options_.kcontrast); - break; - case CHARBONNIER: - charbonnier_diffusivity(evolution_[i].Lx, evolution_[i].Ly, evolution_[i].Lflow, options_.kcontrast); - break; - default: - cerr << "Diffusivity: " << options_.diffusivity << " is not supported" << endl; + // First compute the kcontrast factor + options_.kcontrast = compute_k_percentile(img, options_.kcontrast_percentile, + 1.0, options_.kcontrast_nbins, 0, 0); + + t2 = cv::getTickCount(); + timing_.kcontrast = 1000.0*(t2 - t1) / cv::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); + options_.kcontrast = options_.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 (options_.diffusivity) { + case PM_G1: + pm_g1(evolution_[i].Lx, evolution_[i].Ly, evolution_[i].Lflow, options_.kcontrast); + break; + case PM_G2: + pm_g2(evolution_[i].Lx, evolution_[i].Ly, evolution_[i].Lflow, options_.kcontrast); + break; + case WEICKERT: + weickert_diffusivity(evolution_[i].Lx, evolution_[i].Ly, evolution_[i].Lflow, options_.kcontrast); + break; + case CHARBONNIER: + charbonnier_diffusivity(evolution_[i].Lx, evolution_[i].Ly, evolution_[i].Lflow, options_.kcontrast); + break; + default: + cerr << "Diffusivity: " << options_.diffusivity << " is not supported" << 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]); + } } - // 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 = cv::getTickCount(); + timing_.scale = 1000.0*(t2 - t1) / cv::getTickFrequency(); - t2 = cv::getTickCount(); - timing_.scale = 1000.0*(t2-t1) / cv::getTickFrequency(); - - return 0; + return 0; } /* ************************************************************************* */ /** * @brief This method selects interesting keypoints through the nonlinear scale space * @param kpts Vector of detected keypoints -*/ + */ void AKAZEFeatures::Feature_Detection(std::vector& kpts) { - double t1 = 0.0, t2 = 0.0; + double t1 = 0.0, t2 = 0.0; - t1 = cv::getTickCount(); + t1 = cv::getTickCount(); - vector().swap(kpts); - Compute_Determinant_Hessian_Response(); - Find_Scale_Space_Extrema(kpts); - Do_Subpixel_Refinement(kpts); + kpts.clear(); - t2 = cv::getTickCount(); - timing_.detector = 1000.0*(t2-t1) / cv::getTickFrequency(); + Compute_Determinant_Hessian_Response(); + Find_Scale_Space_Extrema(kpts); + Do_Subpixel_Refinement(kpts); + + t2 = cv::getTickCount(); + timing_.detector = 1000.0*(t2 - t1) / cv::getTickFrequency(); } /* ************************************************************************* */ /** * @brief This method computes the multiscale derivatives for the nonlinear scale space -*/ + */ void AKAZEFeatures::Compute_Multiscale_Derivatives(void) { - double t1 = 0.0, t2 = 0.0; + double t1 = 0.0, t2 = 0.0; - t1 = cv::getTickCount(); + t1 = cv::getTickCount(); #ifdef _OPENMP #pragma omp parallel for #endif - for (int i = 0; i < (int)(evolution_.size()); i++) { + for (int i = 0; i < (int)(evolution_.size()); i++) { - float ratio = pow(2.f,(float)evolution_[i].octave); - int sigma_size_ = fRound(evolution_[i].esigma*options_.derivative_factor/ratio); + float ratio = pow(2.f, (float)evolution_[i].octave); + int sigma_size_ = fRound(evolution_[i].esigma*options_.derivative_factor / 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_); + 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_)); - } + 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 = cv::getTickCount(); - timing_.derivatives = 1000.0*(t2-t1) / cv::getTickFrequency(); + t2 = cv::getTickCount(); + timing_.derivatives = 1000.0*(t2 - t1) / cv::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 AKAZEFeatures::Compute_Determinant_Hessian_Response(void) { - // Firstly compute the multiscale derivatives - Compute_Multiscale_Derivatives(); + // Firstly compute the multiscale derivatives + Compute_Multiscale_Derivatives(); - for (size_t i = 0; i < evolution_.size(); i++) { - if (options_.verbosity == true) { - cout << "Computing detector response. Determinant of Hessian. Evolution time: " << evolution_[i].etime << endl; - } + for (size_t i = 0; i < evolution_.size(); i++) { + if (options_.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); - } + 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 AKAZEFeatures::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; - vector kpts_aux; + 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; + vector kpts_aux; - // Set maximum size - if (options_.descriptor == SURF_UPRIGHT || options_.descriptor == SURF || - options_.descriptor == MLDB_UPRIGHT || options_.descriptor == MLDB) { - smax = 10.0*sqrtf(2.0); - } - else if (options_.descriptor == MSURF_UPRIGHT || options_.descriptor == MSURF) { - smax = 12.0*sqrtf(2.0); - } - - t1 = cv::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 > options_.dthreshold && value >= options_.min_dthreshold && - 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*options_.derivative_factor; - 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; - - // Compare response with the same and lower scale - for (size_t ik = 0; ik < kpts_aux.size(); ik++) { - - if ((point.class_id-1) == kpts_aux[ik].class_id || - point.class_id == kpts_aux[ik].class_id) { - dist = sqrt(pow(point.pt.x*ratio-kpts_aux[ik].pt.x,2)+pow(point.pt.y*ratio-kpts_aux[ik].pt.y,2)); - if (dist <= point.size) { - if (point.response > kpts_aux[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_aux.push_back(point); - npoints++; - } - else { - point.pt.x *= ratio; - point.pt.y *= ratio; - kpts_aux[id_repeated] = point; - } - } // if is_out - } //if is_extremum - } - } // for jx - } // for ix - } // for i - - // Now filter points with the upper scale level - for (size_t i = 0; i < kpts_aux.size(); i++) { - - is_repeated = false; - const cv::KeyPoint& point = kpts_aux[i]; - for (size_t j = i+1; j < kpts_aux.size(); j++) { - - // Compare response with the upper scale - if ((point.class_id+1) == kpts_aux[j].class_id) { - dist = sqrt(pow(point.pt.x-kpts_aux[j].pt.x,2)+pow(point.pt.y-kpts_aux[j].pt.y,2)); - if (dist <= point.size) { - if (point.response < kpts_aux[j].response) { - is_repeated = true; - break; - } - } - } + // Set maximum size + if (options_.descriptor == SURF_UPRIGHT || options_.descriptor == SURF || + options_.descriptor == MLDB_UPRIGHT || options_.descriptor == MLDB) { + smax = 10.0*sqrtf(2.0); + } + else if (options_.descriptor == MSURF_UPRIGHT || options_.descriptor == MSURF) { + smax = 12.0*sqrtf(2.0); } - if (is_repeated == false) - kpts.push_back(point); - } + t1 = cv::getTickCount(); - t2 = cv::getTickCount(); - timing_.extrema = 1000.0*(t2-t1) / cv::getTickFrequency(); + 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 > options_.dthreshold && value >= options_.min_dthreshold && + 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*options_.derivative_factor; + 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; + + // Compare response with the same and lower scale + for (size_t ik = 0; ik < kpts_aux.size(); ik++) { + + if ((point.class_id - 1) == kpts_aux[ik].class_id || + point.class_id == kpts_aux[ik].class_id) { + dist = sqrt(pow(point.pt.x*ratio - kpts_aux[ik].pt.x, 2) + pow(point.pt.y*ratio - kpts_aux[ik].pt.y, 2)); + if (dist <= point.size) { + if (point.response > kpts_aux[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_aux.push_back(point); + npoints++; + } + else { + point.pt.x *= ratio; + point.pt.y *= ratio; + kpts_aux[id_repeated] = point; + } + } // if is_out + } //if is_extremum + } + } // for jx + } // for ix + } // for i + + // Now filter points with the upper scale level + for (size_t i = 0; i < kpts_aux.size(); i++) { + + is_repeated = false; + const cv::KeyPoint& point = kpts_aux[i]; + for (size_t j = i + 1; j < kpts_aux.size(); j++) { + + // Compare response with the upper scale + if ((point.class_id + 1) == kpts_aux[j].class_id) { + dist = sqrt(pow(point.pt.x - kpts_aux[j].pt.x, 2) + pow(point.pt.y - kpts_aux[j].pt.y, 2)); + if (dist <= point.size) { + if (point.response < kpts_aux[j].response) { + is_repeated = true; + break; + } + } + } + } + + if (is_repeated == false) + kpts.push_back(point); + } + + t2 = cv::getTickCount(); + timing_.extrema = 1000.0*(t2 - t1) / cv::getTickFrequency(); } /* ************************************************************************* */ /** * @brief This method performs subpixel refinement of the detected keypoints * @param kpts Vector of detected keypoints -*/ + */ void AKAZEFeatures::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; - cv::Mat A = cv::Mat::zeros(2, 2, CV_32F); - cv::Mat b = cv::Mat::zeros(2, 1, CV_32F); - cv::Mat dst = cv::Mat::zeros(2, 1, CV_32F); + 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; + cv::Mat A = cv::Mat::zeros(2, 2, CV_32F); + cv::Mat b = cv::Mat::zeros(2, 1, CV_32F); + cv::Mat dst = cv::Mat::zeros(2, 1, CV_32F); - t1 = cv::getTickCount(); + t1 = cv::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); + 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 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))); + // 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))); + 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))); + 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 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; - cv::solve(A, b, dst, DECOMP_LU); + cv::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; + 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; + // 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--; + } } - // Delete the point since its not stable - else { - kpts.erase(kpts.begin()+i); - i--; - } - } - t2 = cv::getTickCount(); - timing_.subpixel = 1000.0*(t2-t1) / cv::getTickFrequency(); + t2 = cv::getTickCount(); + timing_.subpixel = 1000.0*(t2 - t1) / cv::getTickFrequency(); } /* ************************************************************************* */ @@ -465,49 +466,49 @@ void AKAZEFeatures::Do_Subpixel_Refinement(std::vector& kpts) { * @brief This method performs feature suppression based on 2D distance * @param kpts Vector of keypoints * @param mdist Maximum distance in pixels -*/ + */ void AKAZEFeatures::Feature_Suppression_Distance(std::vector& kpts, float mdist) const { - vector aux; - vector to_delete; - float dist = 0.0, x1 = 0.0, y1 = 0.0, x2 = 0.0, y2 = 0.0; - bool found = false; + 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); + 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; + } + } } - 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]); } - } } - } - 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(); + kpts.clear(); + kpts = aux; + aux.clear(); } /* ************************************************************************* */ @@ -515,104 +516,104 @@ void AKAZEFeatures::Feature_Suppression_Distance(std::vector& kpts * @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 AKAZEFeatures::Compute_Descriptors(std::vector& kpts, cv::Mat& desc) { - double t1 = 0.0, t2 = 0.0; + double t1 = 0.0, t2 = 0.0; - t1 = cv::getTickCount(); + t1 = cv::getTickCount(); - // Allocate memory for the matrix with the descriptors - if (options_.descriptor < MLDB_UPRIGHT) { - desc = cv::Mat::zeros(kpts.size(), 64, CV_32FC1); - } - else { - // We use the full length binary descriptor -> 486 bits - if (options_.descriptor_size == 0) { - int t = (6+36+120)*options_.descriptor_channels; - desc = cv::Mat::zeros(kpts.size(), ceil(t/8.), CV_8UC1); + // Allocate memory for the matrix with the descriptors + if (options_.descriptor < MLDB_UPRIGHT) { + desc = cv::Mat::zeros(kpts.size(), 64, CV_32FC1); } else { - // We use the random bit selection length binary descriptor - desc = cv::Mat::zeros(kpts.size(), ceil(options_.descriptor_size/8.), CV_8UC1); + // We use the full length binary descriptor -> 486 bits + if (options_.descriptor_size == 0) { + int t = (6 + 36 + 120)*options_.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(options_.descriptor_size / 8.), CV_8UC1); + } } - } - switch (options_.descriptor) { + switch (options_.descriptor) { - case SURF_UPRIGHT : // Upright descriptors, not invariant to rotation + case SURF_UPRIGHT: // Upright descriptors, not invariant to rotation { #ifdef _OPENMP #pragma omp parallel for #endif - for (int i = 0; i < (int)(kpts.size()); i++) { - Get_SURF_Descriptor_Upright_64(kpts[i],desc.ptr(i)); - } + for (int i = 0; i < (int)(kpts.size()); i++) { + Get_SURF_Descriptor_Upright_64(kpts[i], desc.ptr(i)); + } } - break; - case SURF : + break; + case SURF: { #ifdef _OPENMP #pragma omp parallel for #endif - for (int i = 0; i < (int)(kpts.size()); i++) { - Compute_Main_Orientation(kpts[i]); - Get_SURF_Descriptor_64(kpts[i],desc.ptr(i)); - } + for (int i = 0; i < (int)(kpts.size()); i++) { + Compute_Main_Orientation(kpts[i]); + Get_SURF_Descriptor_64(kpts[i], desc.ptr(i)); + } } - break; - case MSURF_UPRIGHT : // Upright descriptors, not invariant to rotation + break; + case MSURF_UPRIGHT: // Upright descriptors, not invariant to rotation { #ifdef _OPENMP #pragma omp parallel for #endif - for (int i = 0; i < (int)(kpts.size()); i++) { - Get_MSURF_Upright_Descriptor_64(kpts[i],desc.ptr(i)); - } + for (int i = 0; i < (int)(kpts.size()); i++) { + Get_MSURF_Upright_Descriptor_64(kpts[i], desc.ptr(i)); + } } - break; - case MSURF : + break; + case MSURF: { #ifdef _OPENMP #pragma omp parallel for #endif - for (int i = 0; i < (int)(kpts.size()); i++) { - Compute_Main_Orientation(kpts[i]); - Get_MSURF_Descriptor_64(kpts[i],desc.ptr(i)); - } + for (int i = 0; i < (int)(kpts.size()); i++) { + Compute_Main_Orientation(kpts[i]); + Get_MSURF_Descriptor_64(kpts[i], desc.ptr(i)); + } } - break; - case MLDB_UPRIGHT : // Upright descriptors, not invariant to rotation + break; + case MLDB_UPRIGHT: // Upright descriptors, not invariant to rotation { #ifdef _OPENMP #pragma omp parallel for #endif - for (int i = 0; i < (int)(kpts.size()); i++) { - if (options_.descriptor_size == 0) - Get_Upright_MLDB_Full_Descriptor(kpts[i], desc.ptr(i)); - else - Get_Upright_MLDB_Descriptor_Subset(kpts[i], desc.ptr(i)); - } + for (int i = 0; i < (int)(kpts.size()); i++) { + if (options_.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 : + break; + case MLDB: { #ifdef _OPENMP #pragma omp parallel for #endif - for (int i = 0; i < (int)(kpts.size()); i++) { - Compute_Main_Orientation(kpts[i]); - if (options_.descriptor_size == 0) - Get_MLDB_Full_Descriptor(kpts[i], desc.ptr(i)); - else - Get_MLDB_Descriptor_Subset(kpts[i], desc.ptr(i)); - } + for (int i = 0; i < (int)(kpts.size()); i++) { + Compute_Main_Orientation(kpts[i]); + if (options_.descriptor_size == 0) + Get_MLDB_Full_Descriptor(kpts[i], desc.ptr(i)); + else + Get_MLDB_Descriptor_Subset(kpts[i], desc.ptr(i)); + } + } + break; } - break; - } - t2 = cv::getTickCount(); - timing_.descriptor = 1000.0*(t2-t1) / cv::getTickFrequency(); + t2 = cv::getTickCount(); + timing_.descriptor = 1000.0*(t2 - t1) / cv::getTickFrequency(); } /* ************************************************************************* */ @@ -621,70 +622,70 @@ void AKAZEFeatures::Compute_Descriptors(std::vector& kpts, cv::Mat * @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 AKAZEFeatures::Compute_Main_Orientation(cv::KeyPoint& kpt) const { - 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}; + 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; + // 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)); + gweight = gauss25[id[i + 6]][id[j + 6]]; + resX[idx] = gweight*(*(evolution_[level].Lx.ptr(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]; - } + Ang[idx] = get_angle(resX[idx], resY[idx]); + ++idx; + } + } } - // 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); + // 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); + } } - } } /* ************************************************************************* */ @@ -694,85 +695,85 @@ void AKAZEFeatures::Compute_Main_Orientation(cv::KeyPoint& kpt) const { * @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 AKAZEFeatures::Get_SURF_Descriptor_Upright_64(const cv::KeyPoint& kpt, float *desc) const { - 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; + 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; + // 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].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; + 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); + // 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; } - } - - // 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); + // convert to unit vector + len = sqrt(len); - for (int i = 0; i < dsize; i++) { - desc[i] /= len; - } + for (int i = 0; i < dsize; i++) { + desc[i] /= len; + } } /* ************************************************************************* */ @@ -784,92 +785,92 @@ void AKAZEFeatures::Get_SURF_Descriptor_Upright_64(const cv::KeyPoint& kpt, floa * @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 AKAZEFeatures::Get_SURF_Descriptor_64(const cv::KeyPoint& kpt, float *desc) const { - 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; + 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; + // 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].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; + 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; + // 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); + // 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; } - } - - // 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); + // convert to unit vector + len = sqrt(len); - for (int i = 0; i < dsize; i++) { - desc[i] /= len; - } + for (int i = 0; i < dsize; i++) { + desc[i] /= len; + } } /* ************************************************************************* */ @@ -881,116 +882,116 @@ void AKAZEFeatures::Get_SURF_Descriptor_64(const cv::KeyPoint& kpt, float *desc) * @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 AKAZEFeatures::Get_MSURF_Upright_Descriptor_64(const cv::KeyPoint& kpt, float *desc) const { - 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; + 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; + // 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; + // 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].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; + 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; + 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); + // 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; } - } - // 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; } - i += 9; - } + // convert to unit vector + len = sqrt(len); - // convert to unit vector - len = sqrt(len); - - for (int i = 0; i < dsize; i++) { - desc[i] /= len; - } + for (int i = 0; i < dsize; i++) { + desc[i] /= len; + } } /* ************************************************************************* */ @@ -1002,120 +1003,120 @@ void AKAZEFeatures::Get_MSURF_Upright_Descriptor_64(const cv::KeyPoint& kpt, flo * @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 AKAZEFeatures::Get_MSURF_Descriptor_64(const cv::KeyPoint& kpt, float *desc) const { - 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; + 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; + // 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; + // 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].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; + 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); + // 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); + // 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; } - } - // 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; } - i += 9; - } + // convert to unit vector + len = sqrt(len); - // convert to unit vector - len = sqrt(len); - - for (int i = 0; i < dsize; i++) { - desc[i] /= len; - } + for (int i = 0; i < dsize; i++) { + desc[i] /= len; + } } /* ************************************************************************* */ @@ -1124,212 +1125,212 @@ void AKAZEFeatures::Get_MSURF_Descriptor_64(const cv::KeyPoint& kpt, float *desc * the provided keypoint * @param kpt Input keypoint * @param desc Descriptor vector -*/ + */ void AKAZEFeatures::Get_Upright_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned char *desc) const { - 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; + 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 - cv::Mat values_1 = cv::Mat::zeros(4, options_.descriptor_channels, CV_32FC1); - cv::Mat values_2 = cv::Mat::zeros(9, options_.descriptor_channels, CV_32FC1); - cv::Mat values_3 = cv::Mat::zeros(16, options_.descriptor_channels, CV_32FC1); + // Matrices for the M-LDB descriptor + cv::Mat values_1 = cv::Mat::zeros(4, options_.descriptor_channels, CV_32FC1); + cv::Mat values_2 = cv::Mat::zeros(9, options_.descriptor_channels, CV_32FC1); + cv::Mat values_3 = cv::Mat::zeros(16, options_.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); + 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 += 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++; } - } - - 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++; + // 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)+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++; + if (*(values_1.ptr(i)+2) > *(values_1.ptr(j)+2)) { + desc[dcount1 / 8] |= (1 << (dcount1 % 8)); + } + dcount1++; } - } - - di /= nsamples; - dx /= nsamples; - dy /= nsamples; - - *(values_2.ptr(dcount2)) = di; - *(values_2.ptr(dcount2)+1) = dx; - *(values_2.ptr(dcount2)+2) = dy; - dcount2++; } - } - //Do binary comparison second level - dcount2 = 0; - 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++; + // Second 3x3 grid + sample_step = ceil(pattern_size*2. / 3.); + dcount2 = 0; - if (*(values_2.ptr(i)+1) > *(values_2.ptr(j)+1)) { - desc[dcount1/8] |= (1<<(dcount1%8)); - } - dcount1++; + 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; - if(*(values_2.ptr(i)+2) > *(values_2.ptr(j)+2)) { - desc[dcount1/8] |= (1<<(dcount1%8)); - } - dcount1++; - } - } + for (int k = i; k < i + sample_step; k++) { + for (int l = j; l < j + sample_step; l++) { - // Third 4x4 grid - sample_step = pattern_size/2; - dcount2 = 0; + // Get the coordinates of the sample point + sample_y = yf + l*scale; + sample_x = xf + k*scale; - 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; + y1 = fRound(sample_y); + x1 = fRound(sample_x); - for (int k = i; k < i + sample_step; k++) { - for (int l = j; l < j + sample_step; l++) { + ri = *(evolution_[level].Lt.ptr(y1)+x1); + rx = *(evolution_[level].Lx.ptr(y1)+x1); + ry = *(evolution_[level].Ly.ptr(y1)+x1); - // Get the coordinates of the sample point - sample_y = yf + l*scale; - sample_x = xf + k*scale; + di += ri; + dx += rx; + dy += ry; + nsamples++; + } + } - y1 = fRound(sample_y); - x1 = fRound(sample_x); + di /= nsamples; + dx /= nsamples; + dy /= nsamples; - 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++; + *(values_2.ptr(dcount2)) = di; + *(values_2.ptr(dcount2)+1) = dx; + *(values_2.ptr(dcount2)+2) = dy; + dcount2++; } - } - - di /= nsamples; - dx /= nsamples; - dy /= nsamples; - - *(values_3.ptr(dcount2)) = di; - *(values_3.ptr(dcount2)+1) = dx; - *(values_3.ptr(dcount2)+2) = dy; - dcount2++; } - } - //Do binary comparison third level - dcount2 = 0; - 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++; + //Do binary comparison second level + dcount2 = 0; + 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_3.ptr(i)+1) > *(values_3.ptr(j)+1)) { - 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_3.ptr(i)+2) > *(values_3.ptr(j)+2)) { - 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++; + } + } + + //Do binary comparison third level + dcount2 = 0; + 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++; + } } - } } /* ************************************************************************* */ @@ -1338,296 +1339,296 @@ void AKAZEFeatures::Get_Upright_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, un * main orientation of the keypoint * @param kpt Input keypoint * @param desc Descriptor vector -*/ + */ void AKAZEFeatures::Get_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned char *desc) const { - 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; + 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 - cv::Mat values_1 = cv::Mat::zeros(4, options_.descriptor_channels, CV_32FC1); - cv::Mat values_2 = cv::Mat::zeros(9, options_.descriptor_channels, CV_32FC1); - cv::Mat values_3 = cv::Mat::zeros(16, options_.descriptor_channels, CV_32FC1); + // Matrices for the M-LDB descriptor + cv::Mat values_1 = cv::Mat::zeros(4, options_.descriptor_channels, CV_32FC1); + cv::Mat values_2 = cv::Mat::zeros(9, options_.descriptor_channels, CV_32FC1); + cv::Mat values_3 = cv::Mat::zeros(16, options_.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); + ri = *(evolution_[level].Lt.ptr(y1)+x1); + rx = *(evolution_[level].Lx.ptr(y1)+x1); + ry = *(evolution_[level].Ly.ptr(y1)+x1); - di += ri; + di += ri; - if (options_.descriptor_channels == 2) { - dx += sqrtf(rx*rx + ry*ry); - } - else if (options_.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; - } + if (options_.descriptor_channels == 2) { + dx += sqrtf(rx*rx + ry*ry); + } + else if (options_.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++; + nsamples++; + } + } + + di /= nsamples; + dx /= nsamples; + dy /= nsamples; + + *(values_1.ptr(dcount2)) = di; + if (options_.descriptor_channels > 1) { + *(values_1.ptr(dcount2)+1) = dx; + } + + if (options_.descriptor_channels > 2) { + *(values_1.ptr(dcount2)+2) = dy; + } + + dcount2++; } - } - - di /= nsamples; - dx /= nsamples; - dy /= nsamples; - - *(values_1.ptr(dcount2)) = di; - if (options_.descriptor_channels > 1 ) { - *(values_1.ptr(dcount2)+1) = dx; - } - - if (options_.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 (options_.descriptor_channels > 1) { + // 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)+1) > *(values_1.ptr(j)+1)) { - desc[dcount1/8] |= (1<<(dcount1%8)); + for (int j = i + 1; j < 4; j++) { + if (*(values_1.ptr(i)) > *(values_1.ptr(j))) { + desc[dcount1 / 8] |= (1 << (dcount1 % 8)); + } + dcount1++; } - - dcount1++; - } } - } - if (options_.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)); + if (options_.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++; + } } - 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 (options_.descriptor_channels == 2) { - dx += sqrtf(rx*rx + ry*ry); - } - else if (options_.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++; + if (options_.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++; + } } - } - - di /= nsamples; - dx /= nsamples; - dy /= nsamples; - - *(values_2.ptr(dcount2)) = di; - if (options_.descriptor_channels > 1) { - *(values_2.ptr(dcount2)+1) = dx; - } - - if (options_.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++; + // 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 (options_.descriptor_channels == 2) { + dx += sqrtf(rx*rx + ry*ry); + } + else if (options_.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 (options_.descriptor_channels > 1) { + *(values_2.ptr(dcount2)+1) = dx; + } + + if (options_.descriptor_channels > 2) { + *(values_2.ptr(dcount2)+2) = dy; + } + + dcount2++; + } } - } - if (options_.descriptor_channels > 1) { + // 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)+1) > *(values_2.ptr(j)+1)) { - desc[dcount1/8] |= (1<<(dcount1%8)); + for (int j = i + 1; j < 9; j++) { + if (*(values_2.ptr(i)) > *(values_2.ptr(j))) { + desc[dcount1 / 8] |= (1 << (dcount1 % 8)); + } + dcount1++; } - dcount1++; - } } - } - if (options_.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)); + if (options_.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++; + } } - 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 (options_.descriptor_channels == 2) { - dx += sqrtf(rx*rx + ry*ry); - } - else if (options_.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++; + if (options_.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++; + } } - } - - di /= nsamples; - dx /= nsamples; - dy /= nsamples; - - *(values_3.ptr(dcount2)) = di; - if (options_.descriptor_channels > 1) - *(values_3.ptr(dcount2)+1) = dx; - - if (options_.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++; + // 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 (options_.descriptor_channels == 2) { + dx += sqrtf(rx*rx + ry*ry); + } + else if (options_.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 (options_.descriptor_channels > 1) + *(values_3.ptr(dcount2)+1) = dx; + + if (options_.descriptor_channels > 2) + *(values_3.ptr(dcount2)+2) = dy; + + dcount2++; + } } - } - if (options_.descriptor_channels > 1) { + // 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)+1) > *(values_3.ptr(j)+1)) { - desc[dcount1/8] |= (1<<(dcount1%8)); + for (int j = i + 1; j < 16; j++) { + if (*(values_3.ptr(i)) > *(values_3.ptr(j))) { + desc[dcount1 / 8] |= (1 << (dcount1 % 8)); + } + dcount1++; } - dcount1++; - } } - } - if (options_.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)); + if (options_.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 (options_.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++; + } } - dcount1++; - } } - } } /* ************************************************************************* */ @@ -1637,88 +1638,88 @@ void AKAZEFeatures::Get_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned c * the bits of the whole descriptor * @param kpt Input keypoint * @param desc Descriptor vector -*/ + */ void AKAZEFeatures::Get_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char *desc) { - float di = 0.f, dx = 0.f, dy = 0.f; - float rx = 0.f, ry = 0.f; - float sample_x = 0.f, sample_y = 0.f; - int x1 = 0, y1 = 0; + float di = 0.f, dx = 0.f, dy = 0.f; + float rx = 0.f, ry = 0.f; + 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)*options_.descriptor_channels, 1); + // Allocate memory for the matrix of values + cv::Mat values = cv::Mat_::zeros((4 + 9 + 16)*options_.descriptor_channels, 1); - // Sample everything, but only do the comparisons - vector steps(3); - steps.at(0) = options_.descriptor_pattern_size; - steps.at(1) = ceil(2.f*options_.descriptor_pattern_size/3.f); - steps.at(2) = options_.descriptor_pattern_size/2; + // Sample everything, but only do the comparisons + vector steps(3); + steps.at(0) = options_.descriptor_pattern_size; + steps.at(1) = ceil(2.f*options_.descriptor_pattern_size / 3.f); + steps.at(2) = options_.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 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++) { + 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); + // 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); + y1 = fRound(sample_y); + x1 = fRound(sample_x); - di += *(evolution_[level].Lt.ptr(y1)+x1); + di += *(evolution_[level].Lt.ptr(y1)+x1); - if (options_.descriptor_channels > 1) { - rx = *(evolution_[level].Lx.ptr(y1)+x1); - ry = *(evolution_[level].Ly.ptr(y1)+x1); + if (options_.descriptor_channels > 1) { + rx = *(evolution_[level].Lx.ptr(y1)+x1); + ry = *(evolution_[level].Ly.ptr(y1)+x1); - if (options_.descriptor_channels == 2) { - dx += sqrtf(rx*rx + ry*ry); - } - else if (options_.descriptor_channels == 3) { - // Get the x and y derivatives on the rotated axis - dx += rx*co + ry*si; - dy += -rx*si + ry*co; - } + if (options_.descriptor_channels == 2) { + dx += sqrtf(rx*rx + ry*ry); + } + else if (options_.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(options_.descriptor_channels*i)) = di; + + if (options_.descriptor_channels == 2) { + *(values.ptr(options_.descriptor_channels*i + 1)) = dx; + } + else if (options_.descriptor_channels == 3) { + *(values.ptr(options_.descriptor_channels*i + 1)) = dx; + *(values.ptr(options_.descriptor_channels*i + 2)) = dy; } - } } - *(values.ptr(options_.descriptor_channels*i)) = di; + // Do the comparisons + const float *vals = values.ptr(0); + const int *comps = descriptorBits_.ptr(0); - if (options_.descriptor_channels == 2) { - *(values.ptr(options_.descriptor_channels*i+1)) = dx; + for (int i = 0; i vals[comps[2 * i + 1]]) { + desc[i / 8] |= (1 << (i % 8)); + } } - else if (options_.descriptor_channels == 3) { - *(values.ptr(options_.descriptor_channels*i+1)) = dx; - *(values.ptr(options_.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)); - } - } } /* ************************************************************************* */ @@ -1728,80 +1729,80 @@ void AKAZEFeatures::Get_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned * The descriptor is computed based on a subset of the bits of the whole descriptor * @param kpt Input keypoint * @param desc Descriptor vector -*/ + */ void AKAZEFeatures::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; + 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)*options_.descriptor_channels, 1); + // Allocate memory for the matrix of values + Mat values = cv::Mat_::zeros((4 + 9 + 16)*options_.descriptor_channels, 1); - vector steps(3); - steps.at(0) = options_.descriptor_pattern_size; - steps.at(1) = ceil(2.f*options_.descriptor_pattern_size/3.f); - steps.at(2) = options_.descriptor_pattern_size/2; + vector steps(3); + steps.at(0) = options_.descriptor_pattern_size; + steps.at(1) = ceil(2.f*options_.descriptor_pattern_size / 3.f); + steps.at(2) = options_.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 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++) { + 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; + // 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); + y1 = fRound(sample_y); + x1 = fRound(sample_x); + di += *(evolution_[level].Lt.ptr(y1)+x1); - if (options_.descriptor_channels > 1) { - rx = *(evolution_[level].Lx.ptr(y1)+x1); - ry = *(evolution_[level].Ly.ptr(y1)+x1); + if (options_.descriptor_channels > 1) { + rx = *(evolution_[level].Lx.ptr(y1)+x1); + ry = *(evolution_[level].Ly.ptr(y1)+x1); - if (options_.descriptor_channels == 2) { - dx += sqrtf(rx*rx + ry*ry); - } - else if (options_.descriptor_channels == 3) { - dx += rx; - dy += ry; - } + if (options_.descriptor_channels == 2) { + dx += sqrtf(rx*rx + ry*ry); + } + else if (options_.descriptor_channels == 3) { + dx += rx; + dy += ry; + } + } + } + } + + *(values.ptr(options_.descriptor_channels*i)) = di; + + if (options_.descriptor_channels == 2) { + *(values.ptr(options_.descriptor_channels*i + 1)) = dx; + } + else if (options_.descriptor_channels == 3) { + *(values.ptr(options_.descriptor_channels*i + 1)) = dx; + *(values.ptr(options_.descriptor_channels*i + 2)) = dy; } - } } - *(values.ptr(options_.descriptor_channels*i)) = di; + // Do the comparisons + const float *vals = values.ptr(0); + const int *comps = descriptorBits_.ptr(0); - if (options_.descriptor_channels == 2) { - *(values.ptr(options_.descriptor_channels*i+1)) = dx; + for (int i = 0; i vals[comps[2 * i + 1]]) { + desc[i / 8] |= (1 << (i % 8)); + } } - else if (options_.descriptor_channels == 3) { - *(values.ptr(options_.descriptor_channels*i+1)) = dx; - *(values.ptr(options_.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)); - } - } } @@ -1809,15 +1810,15 @@ void AKAZEFeatures::Get_Upright_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, /* ************************************************************************* */ /** * @brief This method displays the computation times -*/ + */ void AKAZEFeatures::Show_Computation_Times() const { - cout << "(*) Time Scale Space: " << timing_.scale << endl; - cout << "(*) Time Detector: " << timing_.detector << endl; - cout << " - Time Derivatives: " << timing_.derivatives << endl; - cout << " - Time Extrema: " << timing_.extrema << endl; - cout << " - Time Subpixel: " << timing_.subpixel << endl; - cout << "(*) Time Descriptor: " << timing_.descriptor << endl; - cout << endl; + cout << "(*) Time Scale Space: " << timing_.scale << endl; + cout << "(*) Time Detector: " << timing_.detector << endl; + cout << " - Time Derivatives: " << timing_.derivatives << endl; + cout << " - Time Extrema: " << timing_.extrema << endl; + cout << " - Time Subpixel: " << timing_.subpixel << endl; + cout << "(*) Time Descriptor: " << timing_.descriptor << endl; + cout << endl; } /* ************************************************************************* */ @@ -1835,131 +1836,131 @@ void AKAZEFeatures::Show_Computation_Times() const { * 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 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; + CV_Assert(nbits <= ssz && "descriptor size can't be bigger than full descriptor"); - // 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; + // 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 - for (size_t i=0; i 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 < gsz; j++) { + for (int k = j + 1; k < gsz; k++, c++) { + fullM(c, 0) = i; + fullM(c, 1) = psz*(j % gdiv) - pattern_size; + fullM(c, 2) = psz*(j / gdiv) - pattern_size; + fullM(c, 3) = psz*(k % gdiv) - pattern_size; + fullM(c, 4) = psz*(k / gdiv) - pattern_size; + } + } } - bool n = true; + srand(1024); + Mat_ comps = Mat_(nchannels*ceil(nbits / (float)nchannels), 2); + comps = 1000; - for (int j=0; j samples(29, 3); + Mat_ fullcopy = fullM.clone(); + samples = -1; + + for (size_t i = 0; i < npicks; i++) { + size_t k = rand() % (fullM.rows - i); + if (i < 6) { + // Force use of the coarser grid values and comparisons + k = i; + } + + bool n = true; + + for (int j = 0; j < count; j++) { + if (samples(j, 0) == fullcopy(k, 0) && samples(j, 1) == fullcopy(k, 1) && samples(j, 2) == fullcopy(k, 2)) { + n = false; + comps(i*nchannels, 0) = nchannels*j; + comps(i*nchannels + 1, 0) = nchannels*j + 1; + comps(i*nchannels + 2, 0) = nchannels*j + 2; + break; + } + } + + if (n) { + samples(count, 0) = fullcopy(k, 0); + samples(count, 1) = fullcopy(k, 1); + samples(count, 2) = fullcopy(k, 2); + comps(i*nchannels, 0) = nchannels*count; + comps(i*nchannels + 1, 0) = nchannels*count + 1; + comps(i*nchannels + 2, 0) = nchannels*count + 2; + count++; + } + + n = true; + for (int j = 0; j < count; j++) { + if (samples(j, 0) == fullcopy(k, 0) && samples(j, 1) == fullcopy(k, 3) && samples(j, 2) == fullcopy(k, 4)) { + n = false; + comps(i*nchannels, 1) = nchannels*j; + comps(i*nchannels + 1, 1) = nchannels*j + 1; + comps(i*nchannels + 2, 1) = nchannels*j + 2; + break; + } + } + + if (n) { + samples(count, 0) = fullcopy(k, 0); + samples(count, 1) = fullcopy(k, 3); + samples(count, 2) = fullcopy(k, 4); + comps(i*nchannels, 1) = nchannels*count; + comps(i*nchannels + 1, 1) = nchannels*count + 1; + comps(i*nchannels + 2, 1) = nchannels*count + 2; + count++; + } + + Mat tmp = fullcopy.row(k); + fullcopy.row(fullcopy.rows - i - 1).copyTo(tmp); } - if (n) { - samples(count,0) = fullcopy(k,0); - samples(count,1) = fullcopy(k,1); - samples(count,2) = fullcopy(k,2); - comps(i*nchannels,0) = nchannels*count; - comps(i*nchannels+1,0) = nchannels*count+1; - comps(i*nchannels+2,0) = nchannels*count+2; - count++; - } - - n = true; - for (int j=0; j= 0 && y >= 0) { - return atanf(y/x); - } + if (x >= 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 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); - } + if (x >= 0 && y < 0) { + return 2.0*CV_PI - atanf(-y / x); + } - return 0; + return 0; } /* ************************************************************************* */ @@ -1968,9 +1969,9 @@ inline float get_angle(float x, float y) { * @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)); + return expf(-(x*x + y*y) / (2.0f*sigma*sigma)); } /* ************************************************************************* */ @@ -1980,24 +1981,24 @@ inline float gaussian(float x, float y, float sigma) { * @param y Y Position * @param width Image width * @param height Image height -*/ + */ inline void check_descriptor_limits(int &x, int &y, int width, int height) { - if (x < 0) { - x = 0; - } + if (x < 0) { + x = 0; + } - if (y < 0) { - y = 0; - } + if (y < 0) { + y = 0; + } - if (x > width-1) { - x = width-1; - } + if (x > width - 1) { + x = width - 1; + } - if (y > height-1) { - y = height-1; - } + if (y > height - 1) { + y = height - 1; + } } /* ************************************************************************* */ @@ -2007,6 +2008,6 @@ inline void check_descriptor_limits(int &x, int &y, int width, int height) { * @return dst Nearest integer */ inline int fRound(float flt) { - return (int)(flt+0.5f); + return (int)(flt + 0.5f); }