From: f-morozov Date: Wed, 30 Jul 2014 14:02:08 +0000 (+0400) Subject: AKAZE fixes, tests and tutorial X-Git-Tag: submit/tizen_ivi/20141117.190038~2^2~247^2 X-Git-Url: http://review.tizen.org/git/?a=commitdiff_plain;h=7f829608973148403027e71f26fb704e7717225f;p=profile%2Fivi%2Fopencv.git AKAZE fixes, tests and tutorial --- diff --git a/doc/tutorials/features2d/akaze_matching/akaze_matching.rst b/doc/tutorials/features2d/akaze_matching/akaze_matching.rst new file mode 100644 index 0000000..3fe5df4 --- /dev/null +++ b/doc/tutorials/features2d/akaze_matching/akaze_matching.rst @@ -0,0 +1,161 @@ +.. _akazeMatching: + + +AKAZE local features matching +****************************** + +Introduction +------------------ + +In this tutorial we will learn how to use [AKAZE]_ local features to detect and match keypoints on two images. + +We will find keypoints on a pair of images with given homography matrix, +match them and count the number of inliers (i. e. matches that fit in the given homography). + +You can find expanded version of this example here: https://github.com/pablofdezalc/test_kaze_akaze_opencv + +.. [AKAZE] Fast Explicit Diffusion for Accelerated Features in Nonlinear Scale Spaces. Pablo F. Alcantarilla, Jesús Nuevo and Adrien Bartoli. In British Machine Vision Conference (BMVC), Bristol, UK, September 2013. + +Data +------------------ +We are going to use images 1 and 3 from *Graffity* sequence of Oxford dataset. + +.. image:: images/graf.png + :height: 200pt + :width: 320pt + :alt: Graffity + :align: center + +Homography is given by a 3 by 3 matrix: + +.. code-block:: none + + 7.6285898e-01 -2.9922929e-01 2.2567123e+02 + 3.3443473e-01 1.0143901e+00 -7.6999973e+01 + 3.4663091e-04 -1.4364524e-05 1.0000000e+00 + +You can find the images (*graf1.png*, *graf3.png*) and homography (*H1to3p.xml*) in *opencv/samples/cpp*. + +Source Code +=========== +.. literalinclude:: ../../../../samples/cpp/tutorial_code/features2D/AKAZE_match.cpp + :language: cpp + :linenos: + :tab-width: 4 + +Explanation +=========== + +1. **Load images and homography** + + .. code-block:: cpp + + Mat img1 = imread("graf1.png", IMREAD_GRAYSCALE); + Mat img2 = imread("graf3.png", IMREAD_GRAYSCALE); + + Mat homography; + FileStorage fs("H1to3p.xml", FileStorage::READ); + fs.getFirstTopLevelNode() >> homography; + + We are loading grayscale images here. Homography is stored in the xml created with FileStorage. + +2. **Detect keypoints and compute descriptors using AKAZE** + + .. code-block:: cpp + + vector kpts1, kpts2; + Mat desc1, desc2; + + AKAZE akaze; + akaze(img1, noArray(), kpts1, desc1); + akaze(img2, noArray(), kpts2, desc2); + + We create AKAZE object and use it's *operator()* functionality. Since we don't need the *mask* parameter, *noArray()* is used. + +3. **Use brute-force matcher to find 2-nn matches** + + .. code-block:: cpp + + BFMatcher matcher(NORM_HAMMING); + vector< vector > nn_matches; + matcher.knnMatch(desc1, desc2, nn_matches, 2); + + We use Hamming distance, because AKAZE uses binary descriptor by default. + +4. **Use 2-nn matches to find correct keypoint matches** + + .. code-block:: cpp + + for(size_t i = 0; i < nn_matches.size(); i++) { + DMatch first = nn_matches[i][0]; + float dist1 = nn_matches[i][0].distance; + float dist2 = nn_matches[i][1].distance; + + if(dist1 < nn_match_ratio * dist2) { + matched1.push_back(kpts1[first.queryIdx]); + matched2.push_back(kpts2[first.trainIdx]); + } + } + + If the closest match is *ratio* closer than the second closest one, then the match is correct. + +5. **Check if our matches fit in the homography model** + + .. code-block:: cpp + + for(int i = 0; i < matched1.size(); i++) { + Mat col = Mat::ones(3, 1, CV_64F); + col.at(0) = matched1[i].pt.x; + col.at(1) = matched1[i].pt.y; + + col = homography * col; + col /= col.at(2); + float dist = sqrt( pow(col.at(0) - matched2[i].pt.x, 2) + + pow(col.at(1) - matched2[i].pt.y, 2)); + + if(dist < inlier_threshold) { + int new_i = inliers1.size(); + inliers1.push_back(matched1[i]); + inliers2.push_back(matched2[i]); + good_matches.push_back(DMatch(new_i, new_i, 0)); + } + } + + If the distance from first keypoint's projection to the second keypoint is less than threshold, then it it fits in the homography. + + We create a new set of matches for the inliers, because it is required by the drawing function. + +6. **Output results** + + .. code-block:: cpp + + Mat res; + drawMatches(img1, inliers1, img2, inliers2, good_matches, res); + imwrite("res.png", res); + ... + + Here we save the resulting image and print some statistics. + +Results +======= + +Found matches +-------------- + +.. image:: images/res.png + :height: 200pt + :width: 320pt + :alt: Matches + :align: center + +A-KAZE Matching Results +-------------------------- +Keypoints 1: 2943 + +Keypoints 2: 3511 + +Matches: 447 + +Inliers: 308 + +Inliers Ratio: 0.689038 diff --git a/doc/tutorials/features2d/akaze_matching/images/graf.png b/doc/tutorials/features2d/akaze_matching/images/graf.png new file mode 100644 index 0000000..d9213bc Binary files /dev/null and b/doc/tutorials/features2d/akaze_matching/images/graf.png differ diff --git a/doc/tutorials/features2d/akaze_matching/images/res.png b/doc/tutorials/features2d/akaze_matching/images/res.png new file mode 100644 index 0000000..23fd5bd Binary files /dev/null and b/doc/tutorials/features2d/akaze_matching/images/res.png differ diff --git a/doc/tutorials/features2d/table_of_content_features2d/images/AKAZE_Match_Tutorial_Cover.png b/doc/tutorials/features2d/table_of_content_features2d/images/AKAZE_Match_Tutorial_Cover.png new file mode 100755 index 0000000..fdf2007 Binary files /dev/null and b/doc/tutorials/features2d/table_of_content_features2d/images/AKAZE_Match_Tutorial_Cover.png differ diff --git a/doc/tutorials/features2d/table_of_content_features2d/table_of_content_features2d.rst b/doc/tutorials/features2d/table_of_content_features2d/table_of_content_features2d.rst index f410780..bb79ca3 100644 --- a/doc/tutorials/features2d/table_of_content_features2d/table_of_content_features2d.rst +++ b/doc/tutorials/features2d/table_of_content_features2d/table_of_content_features2d.rst @@ -183,6 +183,25 @@ Learn about how to use the feature points detectors, descriptors and matching f :height: 90pt :width: 90pt ++ + .. tabularcolumns:: m{100pt} m{300pt} + .. cssclass:: toctableopencv + + ===================== ============================================== + |AkazeMatch| **Title:** :ref:`akazeMatching` + + *Compatibility:* > OpenCV 3.0 + + *Author:* Fedor Morozov + + Use *AKAZE* local features to find correspondence between two images. + + ===================== ============================================== + + .. |AkazeMatch| image:: images/AKAZE_Match_Tutorial_Cover.png + :height: 90pt + :width: 90pt + .. raw:: latex \pagebreak @@ -201,3 +220,4 @@ Learn about how to use the feature points detectors, descriptors and matching f ../feature_flann_matcher/feature_flann_matcher ../feature_homography/feature_homography ../detection_of_planar_objects/detection_of_planar_objects + ../akaze_matching/akaze_matching diff --git a/modules/features2d/doc/feature_detection_and_description.rst b/modules/features2d/doc/feature_detection_and_description.rst index 409fe54..4aa7dd3 100644 --- a/modules/features2d/doc/feature_detection_and_description.rst +++ b/modules/features2d/doc/feature_detection_and_description.rst @@ -31,7 +31,7 @@ Detects corners using the FAST algorithm Detects corners using the FAST algorithm by [Rosten06]_. -..note:: In Python API, types are given as ``cv2.FAST_FEATURE_DETECTOR_TYPE_5_8``, ``cv2.FAST_FEATURE_DETECTOR_TYPE_7_12`` and ``cv2.FAST_FEATURE_DETECTOR_TYPE_9_16``. For corner detection, use ``cv2.FAST.detect()`` method. +.. note:: In Python API, types are given as ``cv2.FAST_FEATURE_DETECTOR_TYPE_5_8``, ``cv2.FAST_FEATURE_DETECTOR_TYPE_7_12`` and ``cv2.FAST_FEATURE_DETECTOR_TYPE_9_16``. For corner detection, use ``cv2.FAST.detect()`` method. .. [Rosten06] E. Rosten. Machine Learning for High-speed Corner Detection, 2006. @@ -254,7 +254,17 @@ KAZE ---- .. ocv:class:: KAZE : public Feature2D -Class implementing the KAZE keypoint detector and descriptor extractor, described in [ABD12]_. +Class implementing the KAZE keypoint detector and descriptor extractor, described in [ABD12]_. :: + + class CV_EXPORTS_W KAZE : public Feature2D + { + public: + CV_WRAP KAZE(); + CV_WRAP explicit KAZE(bool extended, bool upright, float threshold = 0.001f, + int octaves = 4, int sublevels = 4, int diffusivity = DIFF_PM_G2); + }; + +.. note:: AKAZE descriptor can only be used with KAZE or AKAZE keypoints .. [ABD12] KAZE Features. Pablo F. Alcantarilla, Adrien Bartoli and Andrew J. Davison. In European Conference on Computer Vision (ECCV), Fiorenze, Italy, October 2012. @@ -262,12 +272,14 @@ KAZE::KAZE ---------- The KAZE constructor -.. ocv:function:: KAZE::KAZE(bool extended, bool upright) +.. ocv:function:: KAZE::KAZE(bool extended, bool upright, float threshold, int octaves, int sublevels, int diffusivity) :param extended: Set to enable extraction of extended (128-byte) descriptor. :param upright: Set to enable use of upright descriptors (non rotation-invariant). - - + :param threshold: Detector response threshold to accept point + :param octaves: Maximum octave evolution of the image + :param sublevels: Default number of sublevels per scale level + :param diffusivity: Diffusivity type. DIFF_PM_G1, DIFF_PM_G2, DIFF_WEICKERT or DIFF_CHARBONNIER AKAZE ----- @@ -278,25 +290,25 @@ Class implementing the AKAZE keypoint detector and descriptor extractor, describ class CV_EXPORTS_W AKAZE : public Feature2D { public: - /// AKAZE Descriptor Type - enum DESCRIPTOR_TYPE { - DESCRIPTOR_KAZE_UPRIGHT = 2, ///< Upright descriptors, not invariant to rotation - DESCRIPTOR_KAZE = 3, - DESCRIPTOR_MLDB_UPRIGHT = 4, ///< Upright descriptors, not invariant to rotation - DESCRIPTOR_MLDB = 5 - }; CV_WRAP AKAZE(); - explicit AKAZE(DESCRIPTOR_TYPE descriptor_type, int descriptor_size = 0, int descriptor_channels = 3); + CV_WRAP explicit AKAZE(int descriptor_type, int descriptor_size = 0, int descriptor_channels = 3, + float threshold = 0.001f, int octaves = 4, int sublevels = 4, int diffusivity = DIFF_PM_G2); }; +.. note:: AKAZE descriptor can only be used with KAZE or AKAZE keypoints + .. [ANB13] Fast Explicit Diffusion for Accelerated Features in Nonlinear Scale Spaces. Pablo F. Alcantarilla, Jesús Nuevo and Adrien Bartoli. In British Machine Vision Conference (BMVC), Bristol, UK, September 2013. AKAZE::AKAZE ------------ The AKAZE constructor -.. ocv:function:: AKAZE::AKAZE(DESCRIPTOR_TYPE descriptor_type, int descriptor_size = 0, int descriptor_channels = 3) +.. ocv:function:: AKAZE::AKAZE(int descriptor_type, int descriptor_size, int descriptor_channels, float threshold, int octaves, int sublevels, int diffusivity) - :param descriptor_type: Type of the extracted descriptor. + :param descriptor_type: Type of the extracted descriptor: DESCRIPTOR_KAZE, DESCRIPTOR_KAZE_UPRIGHT, DESCRIPTOR_MLDB or DESCRIPTOR_MLDB_UPRIGHT. :param descriptor_size: Size of the descriptor in bits. 0 -> Full size - :param descriptor_channels: Number of channels in the descriptor (1, 2, 3). + :param descriptor_channels: Number of channels in the descriptor (1, 2, 3) + :param threshold: Detector response threshold to accept point + :param octaves: Maximum octave evolution of the image + :param sublevels: Default number of sublevels per scale level + :param diffusivity: Diffusivity type. DIFF_PM_G1, DIFF_PM_G2, DIFF_WEICKERT or DIFF_CHARBONNIER diff --git a/modules/features2d/include/opencv2/features2d.hpp b/modules/features2d/include/opencv2/features2d.hpp index 9f46ee2..5aeab1c 100644 --- a/modules/features2d/include/opencv2/features2d.hpp +++ b/modules/features2d/include/opencv2/features2d.hpp @@ -895,6 +895,22 @@ protected: PixelTestFn test_fn_; }; +// KAZE/AKAZE diffusivity +enum { + DIFF_PM_G1 = 0, + DIFF_PM_G2 = 1, + DIFF_WEICKERT = 2, + DIFF_CHARBONNIER = 3 +}; + +// AKAZE descriptor type +enum { + DESCRIPTOR_KAZE_UPRIGHT = 2, ///< Upright descriptors, not invariant to rotation + DESCRIPTOR_KAZE = 3, + DESCRIPTOR_MLDB_UPRIGHT = 4, ///< Upright descriptors, not invariant to rotation + DESCRIPTOR_MLDB = 5 +}; + /*! KAZE implementation */ @@ -902,7 +918,8 @@ class CV_EXPORTS_W KAZE : public Feature2D { public: CV_WRAP KAZE(); - CV_WRAP explicit KAZE(bool extended, bool upright); + CV_WRAP explicit KAZE(bool extended, bool upright, float threshold = 0.001f, + int octaves = 4, int sublevels = 4, int diffusivity = DIFF_PM_G2); virtual ~KAZE(); @@ -928,6 +945,10 @@ protected: CV_PROP bool extended; CV_PROP bool upright; + CV_PROP float threshold; + CV_PROP int octaves; + CV_PROP int sublevels; + CV_PROP int diffusivity; }; /*! @@ -936,16 +957,9 @@ AKAZE implementation class CV_EXPORTS_W AKAZE : public Feature2D { public: - /// AKAZE Descriptor Type - enum DESCRIPTOR_TYPE { - DESCRIPTOR_KAZE_UPRIGHT = 2, ///< Upright descriptors, not invariant to rotation - DESCRIPTOR_KAZE = 3, - DESCRIPTOR_MLDB_UPRIGHT = 4, ///< Upright descriptors, not invariant to rotation - DESCRIPTOR_MLDB = 5 - }; - CV_WRAP AKAZE(); - explicit AKAZE(DESCRIPTOR_TYPE descriptor_type, int descriptor_size = 0, int descriptor_channels = 3); + CV_WRAP explicit AKAZE(int descriptor_type, int descriptor_size = 0, int descriptor_channels = 3, + float threshold = 0.001f, int octaves = 4, int sublevels = 4, int diffusivity = DIFF_PM_G2); virtual ~AKAZE(); @@ -973,7 +987,10 @@ protected: CV_PROP int descriptor; CV_PROP int descriptor_channels; CV_PROP int descriptor_size; - + CV_PROP float threshold; + CV_PROP int octaves; + CV_PROP int sublevels; + CV_PROP int diffusivity; }; /****************************************************************************************\ * Distance * diff --git a/modules/features2d/src/akaze.cpp b/modules/features2d/src/akaze.cpp index 4b1eb19..1d09d06 100644 --- a/modules/features2d/src/akaze.cpp +++ b/modules/features2d/src/akaze.cpp @@ -49,7 +49,10 @@ http://www.robesafe.com/personal/pablo.alcantarilla/papers/Alcantarilla13bmvc.pd */ #include "precomp.hpp" -#include "akaze/AKAZEFeatures.h" +#include "kaze/AKAZEFeatures.h" + +#include +using namespace std; namespace cv { @@ -57,13 +60,22 @@ namespace cv : descriptor(DESCRIPTOR_MLDB) , descriptor_channels(3) , descriptor_size(0) + , threshold(0.001f) + , octaves(4) + , sublevels(4) + , diffusivity(DIFF_PM_G2) { } - AKAZE::AKAZE(DESCRIPTOR_TYPE _descriptor_type, int _descriptor_size, int _descriptor_channels) + AKAZE::AKAZE(int _descriptor_type, int _descriptor_size, int _descriptor_channels, + float _threshold, int _octaves, int _sublevels, int _diffusivity) : descriptor(_descriptor_type) , descriptor_channels(_descriptor_channels) , descriptor_size(_descriptor_size) + , threshold(_threshold) + , octaves(_octaves) + , sublevels(_sublevels) + , diffusivity(_diffusivity) { } @@ -78,12 +90,12 @@ namespace cv { switch (descriptor) { - case cv::AKAZE::DESCRIPTOR_KAZE: - case cv::AKAZE::DESCRIPTOR_KAZE_UPRIGHT: + case cv::DESCRIPTOR_KAZE: + case cv::DESCRIPTOR_KAZE_UPRIGHT: return 64; - case cv::AKAZE::DESCRIPTOR_MLDB: - case cv::AKAZE::DESCRIPTOR_MLDB_UPRIGHT: + case cv::DESCRIPTOR_MLDB: + case cv::DESCRIPTOR_MLDB_UPRIGHT: // We use the full length binary descriptor -> 486 bits if (descriptor_size == 0) { @@ -106,12 +118,12 @@ namespace cv { switch (descriptor) { - case cv::AKAZE::DESCRIPTOR_KAZE: - case cv::AKAZE::DESCRIPTOR_KAZE_UPRIGHT: + case cv::DESCRIPTOR_KAZE: + case cv::DESCRIPTOR_KAZE_UPRIGHT: return CV_32F; - case cv::AKAZE::DESCRIPTOR_MLDB: - case cv::AKAZE::DESCRIPTOR_MLDB_UPRIGHT: + case cv::DESCRIPTOR_MLDB: + case cv::DESCRIPTOR_MLDB_UPRIGHT: return CV_8U; default: @@ -124,12 +136,12 @@ namespace cv { switch (descriptor) { - case cv::AKAZE::DESCRIPTOR_KAZE: - case cv::AKAZE::DESCRIPTOR_KAZE_UPRIGHT: + case cv::DESCRIPTOR_KAZE: + case cv::DESCRIPTOR_KAZE_UPRIGHT: return cv::NORM_L2; - case cv::AKAZE::DESCRIPTOR_MLDB: - case cv::AKAZE::DESCRIPTOR_MLDB_UPRIGHT: + case cv::DESCRIPTOR_MLDB: + case cv::DESCRIPTOR_MLDB_UPRIGHT: return cv::NORM_HAMMING; default: @@ -153,11 +165,15 @@ namespace cv cv::Mat& desc = descriptors.getMatRef(); AKAZEOptions options; - options.descriptor = static_cast(descriptor); + options.descriptor = descriptor; options.descriptor_channels = descriptor_channels; options.descriptor_size = descriptor_size; options.img_width = img.cols; options.img_height = img.rows; + options.dthreshold = threshold; + options.omax = octaves; + options.nsublevels = sublevels; + options.diffusivity = diffusivity; AKAZEFeatures impl(options); impl.Create_Nonlinear_Scale_Space(img1_32); @@ -188,7 +204,7 @@ namespace cv img.convertTo(img1_32, CV_32F, 1.0 / 255.0, 0); AKAZEOptions options; - options.descriptor = static_cast(descriptor); + options.descriptor = descriptor; options.descriptor_channels = descriptor_channels; options.descriptor_size = descriptor_size; options.img_width = img.cols; @@ -216,7 +232,7 @@ namespace cv cv::Mat& desc = descriptors.getMatRef(); AKAZEOptions options; - options.descriptor = static_cast(descriptor); + options.descriptor = descriptor; options.descriptor_channels = descriptor_channels; options.descriptor_size = descriptor_size; options.img_width = img.cols; @@ -229,4 +245,4 @@ namespace cv CV_Assert((!desc.rows || desc.cols == descriptorSize())); CV_Assert((!desc.rows || (desc.type() == descriptorType()))); } -} \ No newline at end of file +} diff --git a/modules/features2d/src/akaze/AKAZEFeatures.cpp b/modules/features2d/src/akaze/AKAZEFeatures.cpp deleted file mode 100644 index e5955b2..0000000 --- a/modules/features2d/src/akaze/AKAZEFeatures.cpp +++ /dev/null @@ -1,1941 +0,0 @@ -/** - * @file AKAZEFeatures.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 "AKAZEFeatures.h" -#include "../kaze/fed.h" -#include "../kaze/nldiffusion_functions.h" - -using namespace std; -using namespace cv; -using namespace cv::details::kaze; - -/* ************************************************************************* */ -/** - * @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; - - if (options_.descriptor_size > 0 && options_.descriptor >= cv::AKAZE::DESCRIPTOR_MLDB_UPRIGHT) - { - generateDescriptorSubsample(descriptorSamples_, descriptorBits_, options_.descriptor_size, - options_.descriptor_pattern_size, options_.descriptor_channels); - } - - Allocate_Memory_Evolution(); -} - -/* ************************************************************************* */ -/** - * @brief This method allocates the memory for the nonlinear diffusion evolution - */ -void AKAZEFeatures::Allocate_Memory_Evolution(void) { - - float rfactor = 0.0f; - 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.0f / 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; - } - - 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.5f*(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.0f; - ttime = evolution_[i].etime - evolution_[i - 1].etime; - naux = fed_tau_by_process_time(ttime, 1, 0.25f, 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 AKAZEFeatures::Create_Nonlinear_Scale_Space(const cv::Mat& img) -{ - CV_Assert(evolution_.size() > 0); - - // 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.0f, options_.kcontrast_nbins, 0, 0); - - // 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.75f; - } - else { - evolution_[i - 1].Lt.copyTo(evolution_[i].Lt); - } - - gaussian_2D_convolution(evolution_[i].Lt, evolution_[i].Lsmooth, 0, 0, 1.0f); - - // 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 AKAZEOptions::PM_G1: - pm_g1(evolution_[i].Lx, evolution_[i].Ly, evolution_[i].Lflow, options_.kcontrast); - break; - case AKAZEOptions::PM_G2: - pm_g2(evolution_[i].Lx, evolution_[i].Ly, evolution_[i].Lflow, options_.kcontrast); - break; - case AKAZEOptions::WEICKERT: - weickert_diffusivity(evolution_[i].Lx, evolution_[i].Ly, evolution_[i].Lflow, options_.kcontrast); - break; - case AKAZEOptions::CHARBONNIER: - charbonnier_diffusivity(evolution_[i].Lx, evolution_[i].Ly, evolution_[i].Lflow, options_.kcontrast); - break; - default: - CV_Error(options_.diffusivity, "Diffusivity is not supported"); - break; - } - - // Perform FED n inner steps - for (int j = 0; j < nsteps_[i - 1]; j++) { - cv::details::kaze::nld_step_scalar(evolution_[i].Lt, evolution_[i].Lflow, evolution_[i].Lstep, tsteps_[i - 1][j]); - } - } - - 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) -{ - kpts.clear(); - - Compute_Determinant_Hessian_Response(); - Find_Scale_Space_Extrema(kpts); - Do_Subpixel_Refinement(kpts); -} - -/* ************************************************************************* */ - -class MultiscaleDerivativesInvoker : public cv::ParallelLoopBody -{ -public: - explicit MultiscaleDerivativesInvoker(std::vector& ev, const AKAZEOptions& opt) - : evolution_(&ev) - , options_(opt) - { - } - - - void operator()(const cv::Range& range) const - { - std::vector& evolution = *evolution_; - - for (int i = range.start; i < range.end; i++) - { - 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_); - - 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_)); - } - } - -private: - std::vector* evolution_; - AKAZEOptions options_; -}; - -/** - * @brief This method computes the multiscale derivatives for the nonlinear scale space - */ -void AKAZEFeatures::Compute_Multiscale_Derivatives(void) -{ - cv::parallel_for_(cv::Range(0, (int)evolution_.size()), - MultiscaleDerivativesInvoker(evolution_, options_)); -} - -/* ************************************************************************* */ -/** - * @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(); - - for (size_t i = 0; i < evolution_.size(); i++) - { - 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) -{ - - 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 == cv::AKAZE::DESCRIPTOR_MLDB_UPRIGHT || options_.descriptor == cv::AKAZE::DESCRIPTOR_MLDB) { - smax = 10.0f*sqrtf(2.0f); - } - else if (options_.descriptor == cv::AKAZE::DESCRIPTOR_KAZE_UPRIGHT || options_.descriptor == cv::AKAZE::DESCRIPTOR_KAZE) { - smax = 12.0f*sqrtf(2.0f); - } - - 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 = (int)evolution_[i].octave; - point.class_id = (int)i; - ratio = pow(2.f, point.octave); - sigma_size_ = fRound(point.size / ratio); - point.pt.x = static_cast(jx); - point.pt.y = static_cast(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 = (int)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& pt = kpts_aux[i]; - for (size_t j = i + 1; j < kpts_aux.size(); j++) { - - // Compare response with the upper scale - if ((pt.class_id + 1) == kpts_aux[j].class_id) { - dist = sqrt(pow(pt.pt.x - kpts_aux[j].pt.x, 2) + pow(pt.pt.y - kpts_aux[j].pt.y, 2)); - if (dist <= pt.size) { - if (pt.response < kpts_aux[j].response) { - is_repeated = true; - break; - } - } - } - } - - if (is_repeated == false) - kpts.push_back(pt); - } -} - -/* ************************************************************************* */ -/** - * @brief This method performs subpixel refinement of the detected keypoints - * @param kpts Vector of detected keypoints - */ -void AKAZEFeatures::Do_Subpixel_Refinement(std::vector& kpts) -{ - 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); - - 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.5f)*(*(evolution_[kpts[i].class_id].Ldet.ptr(y)+x + 1) - - *(evolution_[kpts[i].class_id].Ldet.ptr(y)+x - 1)); - Dy = (0.5f)*(*(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.0f*(*(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.0f*(*(evolution_[kpts[i].class_id].Ldet.ptr(y)+x))); - - Dxy = (0.25f)*(*(evolution_[kpts[i].class_id].Ldet.ptr(y + 1) + x + 1) - + (*(evolution_[kpts[i].class_id].Ldet.ptr(y - 1) + x - 1))) - - (0.25f)*(*(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; - - cv::solve(A, b, dst, DECOMP_LU); - - if (fabs(*(dst.ptr(0))) <= 1.0f && fabs(*(dst.ptr(1))) <= 1.0f) { - kpts[i].pt.x = x + (*(dst.ptr(0))); - kpts[i].pt.y = y + (*(dst.ptr(1))); - kpts[i].pt.x *= powf(2.f, (float)evolution_[kpts[i].class_id].octave); - kpts[i].pt.y *= powf(2.f, (float)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.0f; - } - // Delete the point since its not stable - else { - kpts.erase(kpts.begin() + i); - i--; - } - } -} - -/* ************************************************************************* */ - -class SURF_Descriptor_Upright_64_Invoker : public cv::ParallelLoopBody -{ -public: - SURF_Descriptor_Upright_64_Invoker(std::vector& kpts, cv::Mat& desc, std::vector& evolution) - : keypoints_(&kpts) - , descriptors_(&desc) - , evolution_(&evolution) - { - } - - void operator() (const Range& range) const - { - for (int i = range.start; i < range.end; i++) - { - Get_SURF_Descriptor_Upright_64((*keypoints_)[i], descriptors_->ptr(i)); - } - } - - void Get_SURF_Descriptor_Upright_64(const cv::KeyPoint& kpt, float* desc) const; - -private: - std::vector* keypoints_; - cv::Mat* descriptors_; - std::vector* evolution_; -}; - -class SURF_Descriptor_64_Invoker : public cv::ParallelLoopBody -{ -public: - SURF_Descriptor_64_Invoker(std::vector& kpts, cv::Mat& desc, std::vector& evolution) - : keypoints_(&kpts) - , descriptors_(&desc) - , evolution_(&evolution) - { - } - - void operator()(const Range& range) const - { - for (int i = range.start; i < range.end; i++) - { - AKAZEFeatures::Compute_Main_Orientation((*keypoints_)[i], *evolution_); - Get_SURF_Descriptor_64((*keypoints_)[i], descriptors_->ptr(i)); - } - } - - void Get_SURF_Descriptor_64(const cv::KeyPoint& kpt, float* desc) const; - -private: - std::vector* keypoints_; - cv::Mat* descriptors_; - std::vector* evolution_; -}; - -class MSURF_Upright_Descriptor_64_Invoker : public cv::ParallelLoopBody -{ -public: - MSURF_Upright_Descriptor_64_Invoker(std::vector& kpts, cv::Mat& desc, std::vector& evolution) - : keypoints_(&kpts) - , descriptors_(&desc) - , evolution_(&evolution) - { - } - - void operator()(const Range& range) const - { - for (int i = range.start; i < range.end; i++) - { - Get_MSURF_Upright_Descriptor_64((*keypoints_)[i], descriptors_->ptr(i)); - } - } - - void Get_MSURF_Upright_Descriptor_64(const cv::KeyPoint& kpt, float* desc) const; - -private: - std::vector* keypoints_; - cv::Mat* descriptors_; - std::vector* evolution_; -}; - -class MSURF_Descriptor_64_Invoker : public cv::ParallelLoopBody -{ -public: - MSURF_Descriptor_64_Invoker(std::vector& kpts, cv::Mat& desc, std::vector& evolution) - : keypoints_(&kpts) - , descriptors_(&desc) - , evolution_(&evolution) - { - } - - void operator() (const Range& range) const - { - for (int i = range.start; i < range.end; i++) - { - AKAZEFeatures::Compute_Main_Orientation((*keypoints_)[i], *evolution_); - Get_MSURF_Descriptor_64((*keypoints_)[i], descriptors_->ptr(i)); - } - } - - void Get_MSURF_Descriptor_64(const cv::KeyPoint& kpt, float* desc) const; - -private: - std::vector* keypoints_; - cv::Mat* descriptors_; - std::vector* evolution_; -}; - -class Upright_MLDB_Full_Descriptor_Invoker : public cv::ParallelLoopBody -{ -public: - Upright_MLDB_Full_Descriptor_Invoker(std::vector& kpts, cv::Mat& desc, std::vector& evolution, AKAZEOptions& options) - : keypoints_(&kpts) - , descriptors_(&desc) - , evolution_(&evolution) - , options_(&options) - { - } - - void operator() (const Range& range) const - { - for (int i = range.start; i < range.end; i++) - { - Get_Upright_MLDB_Full_Descriptor((*keypoints_)[i], descriptors_->ptr(i)); - } - } - - void Get_Upright_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned char* desc) const; - -private: - std::vector* keypoints_; - cv::Mat* descriptors_; - std::vector* evolution_; - AKAZEOptions* options_; -}; - -class Upright_MLDB_Descriptor_Subset_Invoker : public cv::ParallelLoopBody -{ -public: - Upright_MLDB_Descriptor_Subset_Invoker(std::vector& kpts, - cv::Mat& desc, - std::vector& evolution, - AKAZEOptions& options, - cv::Mat descriptorSamples, - cv::Mat descriptorBits) - : keypoints_(&kpts) - , descriptors_(&desc) - , evolution_(&evolution) - , options_(&options) - , descriptorSamples_(descriptorSamples) - , descriptorBits_(descriptorBits) - { - } - - void operator() (const Range& range) const - { - for (int i = range.start; i < range.end; i++) - { - Get_Upright_MLDB_Descriptor_Subset((*keypoints_)[i], descriptors_->ptr(i)); - } - } - - void Get_Upright_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char* desc) const; - -private: - std::vector* keypoints_; - cv::Mat* descriptors_; - std::vector* evolution_; - AKAZEOptions* options_; - - cv::Mat descriptorSamples_; // List of positions in the grids to sample LDB bits from. - cv::Mat descriptorBits_; -}; - -class MLDB_Full_Descriptor_Invoker : public cv::ParallelLoopBody -{ -public: - MLDB_Full_Descriptor_Invoker(std::vector& kpts, cv::Mat& desc, std::vector& evolution, AKAZEOptions& options) - : keypoints_(&kpts) - , descriptors_(&desc) - , evolution_(&evolution) - , options_(&options) - { - } - - void operator() (const Range& range) const - { - for (int i = range.start; i < range.end; i++) - { - AKAZEFeatures::Compute_Main_Orientation((*keypoints_)[i], *evolution_); - Get_MLDB_Full_Descriptor((*keypoints_)[i], descriptors_->ptr(i)); - } - } - - void Get_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned char* desc) const; - -private: - std::vector* keypoints_; - cv::Mat* descriptors_; - std::vector* evolution_; - AKAZEOptions* options_; -}; - -class MLDB_Descriptor_Subset_Invoker : public cv::ParallelLoopBody -{ -public: - MLDB_Descriptor_Subset_Invoker(std::vector& kpts, - cv::Mat& desc, - std::vector& evolution, - AKAZEOptions& options, - cv::Mat descriptorSamples, - cv::Mat descriptorBits) - : keypoints_(&kpts) - , descriptors_(&desc) - , evolution_(&evolution) - , options_(&options) - , descriptorSamples_(descriptorSamples) - , descriptorBits_(descriptorBits) - { - } - - void operator() (const Range& range) const - { - for (int i = range.start; i < range.end; i++) - { - AKAZEFeatures::Compute_Main_Orientation((*keypoints_)[i], *evolution_); - Get_MLDB_Descriptor_Subset((*keypoints_)[i], descriptors_->ptr(i)); - } - } - - void Get_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char* desc) const; - -private: - std::vector* keypoints_; - cv::Mat* descriptors_; - std::vector* evolution_; - AKAZEOptions* options_; - - cv::Mat descriptorSamples_; // List of positions in the grids to sample LDB bits from. - cv::Mat descriptorBits_; -}; - -/** - * @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) -{ - // Allocate memory for the matrix with the descriptors - if (options_.descriptor < cv::AKAZE::DESCRIPTOR_MLDB_UPRIGHT) { - desc = cv::Mat::zeros((int)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((int)kpts.size(), (int)ceil(t / 8.), CV_8UC1); - } - else { - // We use the random bit selection length binary descriptor - desc = cv::Mat::zeros((int)kpts.size(), (int)ceil(options_.descriptor_size / 8.), CV_8UC1); - } - } - - switch (options_.descriptor) - { - case cv::AKAZE::DESCRIPTOR_KAZE_UPRIGHT: // Upright descriptors, not invariant to rotation - { - cv::parallel_for_(cv::Range(0, (int)kpts.size()), MSURF_Upright_Descriptor_64_Invoker(kpts, desc, evolution_)); - } - break; - case cv::AKAZE::DESCRIPTOR_KAZE: - { - cv::parallel_for_(cv::Range(0, (int)kpts.size()), MSURF_Descriptor_64_Invoker(kpts, desc, evolution_)); - } - break; - case cv::AKAZE::DESCRIPTOR_MLDB_UPRIGHT: // Upright descriptors, not invariant to rotation - { - if (options_.descriptor_size == 0) - cv::parallel_for_(cv::Range(0, (int)kpts.size()), Upright_MLDB_Full_Descriptor_Invoker(kpts, desc, evolution_, options_)); - else - cv::parallel_for_(cv::Range(0, (int)kpts.size()), Upright_MLDB_Descriptor_Subset_Invoker(kpts, desc, evolution_, options_, descriptorSamples_, descriptorBits_)); - } - break; - case cv::AKAZE::DESCRIPTOR_MLDB: - { - if (options_.descriptor_size == 0) - cv::parallel_for_(cv::Range(0, (int)kpts.size()), MLDB_Full_Descriptor_Invoker(kpts, desc, evolution_, options_)); - else - cv::parallel_for_(cv::Range(0, (int)kpts.size()), MLDB_Descriptor_Subset_Invoker(kpts, desc, evolution_, options_, descriptorSamples_, descriptorBits_)); - } - break; - } -} - -/* ************************************************************************* */ -/** - * @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 AKAZEFeatures::Compute_Main_Orientation(cv::KeyPoint& kpt, const std::vector& evolution_) { - - 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 << evolution_[level].octave); - s = fRound(0.5f*kpt.size / ratio); - xf = kpt.pt.x / ratio; - yf = kpt.pt.y / ratio; - - // 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); - - 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 < (float)(2.0 * CV_PI); ang1 += 0.15f) { - ang2 = (ang1 + (float)(CV_PI / 3.0) >(float)(2.0*CV_PI) ? ang1 - (float)(5.0*CV_PI / 3.0) : ang1 + (float)(CV_PI / 3.0)); - 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.0f*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 (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 MSURF_Upright_Descriptor_64_Invoker::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; - - // Subregion centers for the 4x4 gaussian weighting - float cx = -0.5f, cy = 0.5f; - - const std::vector& evolution = *evolution_; - - // 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 << kpt.octave); - scale = fRound(0.5f*kpt.size / ratio); - level = kpt.class_id; - yf = kpt.pt.y / ratio; - xf = kpt.pt.x / ratio; - - 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.0f; - cy = -0.5f; - - while (j < pattern_size) { - dx = dy = mdx = mdy = 0.0; - cy += 1.0f; - 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.50f*scale); - - y1 = (int)(sample_y - .5); - x1 = (int)(sample_x - .5); - - y2 = (int)(sample_y + .5); - x2 = (int)(sample_x + .5); - - 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.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + (1.0f - 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.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + (1.0f - 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 (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 MSURF_Descriptor_64_Invoker::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; - - // Subregion centers for the 4x4 gaussian weighting - float cx = -0.5f, cy = 0.5f; - - const std::vector& evolution = *evolution_; - - // 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 << kpt.octave); - scale = fRound(0.5f*kpt.size / ratio); - angle = kpt.angle; - level = kpt.class_id; - yf = kpt.pt.y / ratio; - xf = kpt.pt.x / ratio; - 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.0f; - cy = -0.5f; - - while (j < pattern_size) { - dx = dy = mdx = mdy = 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.5f*scale); - - y1 = fRound(sample_y - 0.5f); - x1 = fRound(sample_x - 0.5f); - - y2 = fRound(sample_y + 0.5f); - x2 = fRound(sample_x + 0.5f); - - 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.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + (1.0f - 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.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + (1.0f - 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 (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 Upright_MLDB_Full_Descriptor_Invoker::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; - - const AKAZEOptions & options = *options_; - const std::vector& evolution = *evolution_; - - // 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 << kpt.octave); - scale = fRound(0.5f*kpt.size / ratio); - level = kpt.class_id; - yf = kpt.pt.y / ratio; - xf = kpt.pt.x / ratio; - - // First 2x2 grid - pattern_size = options_->descriptor_pattern_size; - sample_step = pattern_size; - - 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_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 = static_cast(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++; - } - } - - //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_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++; - } - } - - //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++; - } - } -} - -/* ************************************************************************* */ -/** - * @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 MLDB_Full_Descriptor_Invoker::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; - - const AKAZEOptions & options = *options_; - const std::vector& evolution = *evolution_; - - // 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 << kpt.octave); - scale = fRound(0.5f*kpt.size / ratio); - angle = kpt.angle; - level = kpt.class_id; - yf = kpt.pt.y / ratio; - xf = kpt.pt.x / ratio; - co = cos(angle); - si = sin(angle); - - // First 2x2 grid - pattern_size = options.descriptor_pattern_size; - sample_step = pattern_size; - - 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 (float k = (float)i; k < i + sample_step; k++) { - for (float l = (float)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_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) { - 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 (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++; - } - } - } - - // Second 3x3 grid - sample_step = static_cast(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++; - } - } - - // 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 (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++; - } - } - } - - 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++; - } - } - } - - // 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++; - } - } - - // 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 (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++; - } - } - } -} - -/* ************************************************************************* */ -/** - * @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 MLDB_Descriptor_Subset_Invoker::Get_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char *desc) const { - - 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; - - const AKAZEOptions & options = *options_; - const std::vector& evolution = *evolution_; - - // Get the information from the keypoint - float ratio = (float)(1 << kpt.octave); - int scale = fRound(0.5f*kpt.size / ratio); - float angle = kpt.angle; - int level = kpt.class_id; - float yf = kpt.pt.y / ratio; - float xf = kpt.pt.x / ratio; - float co = cos(angle); - float si = sin(angle); - - // 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) = (int)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++) { - const 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*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 (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; - } - } - } - } - - *(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; - } - } - - // 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 Upright_MLDB_Descriptor_Subset_Invoker::Get_Upright_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char *desc) const { - - 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; - - const AKAZEOptions & options = *options_; - const std::vector& evolution = *evolution_; - - // Get the information from the keypoint - float ratio = (float)(1 << kpt.octave); - int scale = fRound(0.5f*kpt.size / ratio); - int level = kpt.class_id; - float yf = kpt.pt.y / ratio; - float xf = kpt.pt.x / ratio; - - // 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) = static_cast(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++) { - const 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 (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; - } - } - } - } - - *(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; - } - } - - // 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 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 (int i = 0, c = 0; i < 3; i++) { - int gdiv = i + 2; //grid divisions, per row - int gsz = gdiv*gdiv; - int psz = (int)ceil(2.f*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; - } - } - } - - srand(1024); - Mat_ comps = Mat_(nchannels * (int)ceil(nbits / (float)nchannels), 2); - comps = 1000; - - // Select some samples. A sample includes all channels - int count = 0; - int npicks = (int)ceil(nbits / (float)nchannels); - Mat_ samples(29, 3); - Mat_ fullcopy = fullM.clone(); - samples = -1; - - for (int i = 0; i < npicks; i++) { - int 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); - } - - sampleList = samples.rowRange(0, count).clone(); - comparisons = comps.rowRange(0, nbits).clone(); -} - -/* ************************************************************************* */ -/** - * @brief This function computes the angle from the vector given by (X Y). From 0 to 2*Pi - */ -inline float get_angle(float x, float y) { - - if (x >= 0 && y >= 0) { - return atanf(y / x); - } - - if (x < 0 && y >= 0) { - return static_cast(CV_PI)-atanf(-y / x); - } - - if (x < 0 && y < 0) { - return static_cast(CV_PI)+atanf(y / x); - } - - if (x >= 0 && y < 0) { - return static_cast(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, int width, 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); -} \ No newline at end of file diff --git a/modules/features2d/src/akaze/AKAZEFeatures.h b/modules/features2d/src/akaze/AKAZEFeatures.h deleted file mode 100644 index 302ef0d..0000000 --- a/modules/features2d/src/akaze/AKAZEFeatures.h +++ /dev/null @@ -1,65 +0,0 @@ -/** - * @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 - */ - -#pragma once - -/* ************************************************************************* */ -// Includes -#include "precomp.hpp" -#include "AKAZEConfig.h" - -/* ************************************************************************* */ -// AKAZE Class Declaration -class AKAZEFeatures { - -private: - - AKAZEOptions options_; ///< Configuration options for AKAZE - 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 - - /// 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_; - -public: - - /// Constructor with input arguments - AKAZEFeatures(const AKAZEOptions& options); - - /// Scale Space methods - void Allocate_Memory_Evolution(); - 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); - - // Feature description methods - void Compute_Descriptors(std::vector& kpts, cv::Mat& desc); - - static void Compute_Main_Orientation(cv::KeyPoint& kpt, const std::vector& evolution_); -}; - -/* ************************************************************************* */ -// Inline functions - -// 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, int width, int height); -int fRound(float flt); diff --git a/modules/features2d/src/kaze.cpp b/modules/features2d/src/kaze.cpp index 88fb999..910ec27 100644 --- a/modules/features2d/src/kaze.cpp +++ b/modules/features2d/src/kaze.cpp @@ -55,12 +55,21 @@ namespace cv KAZE::KAZE() : extended(false) , upright(false) + , threshold(0.001f) + , octaves(4) + , sublevels(4) + , diffusivity(DIFF_PM_G2) { } - KAZE::KAZE(bool _extended, bool _upright) + KAZE::KAZE(bool _extended, bool _upright, float _threshold, int _octaves, + int _sublevels, int _diffusivity) : extended(_extended) , upright(_upright) + , threshold(_threshold) + , octaves(_octaves) + , sublevels(_sublevels) + , diffusivity(_diffusivity) { } @@ -111,6 +120,10 @@ namespace cv options.img_height = img.rows; options.extended = extended; options.upright = upright; + options.dthreshold = threshold; + options.omax = octaves; + options.nsublevels = sublevels; + options.diffusivity = diffusivity; KAZEFeatures impl(options); impl.Create_Nonlinear_Scale_Space(img1_32); @@ -180,4 +193,4 @@ namespace cv CV_Assert((!desc.rows || desc.cols == descriptorSize())); CV_Assert((!desc.rows || (desc.type() == descriptorType()))); } -} \ No newline at end of file +} diff --git a/modules/features2d/src/akaze/AKAZEConfig.h b/modules/features2d/src/kaze/AKAZEConfig.h similarity index 70% rename from modules/features2d/src/akaze/AKAZEConfig.h rename to modules/features2d/src/kaze/AKAZEConfig.h index 1c1203f..c7ac1cf 100644 --- a/modules/features2d/src/akaze/AKAZEConfig.h +++ b/modules/features2d/src/kaze/AKAZEConfig.h @@ -5,7 +5,8 @@ * @author Pablo F. Alcantarilla, Jesus Nuevo */ -#pragma once +#ifndef __OPENCV_FEATURES_2D_AKAZE_CONFIG_H__ +#define __OPENCV_FEATURES_2D_AKAZE_CONFIG_H__ /* ************************************************************************* */ // OpenCV @@ -28,14 +29,6 @@ const float gauss25[7][7] = { /// AKAZE configuration options structure struct AKAZEOptions { - /// AKAZE Diffusivities - enum DIFFUSIVITY_TYPE { - PM_G1 = 0, - PM_G2 = 1, - WEICKERT = 2, - CHARBONNIER = 3 - }; - AKAZEOptions() : omax(4) , nsublevels(4) @@ -44,12 +37,12 @@ struct AKAZEOptions { , soffset(1.6f) , derivative_factor(1.5f) , sderivatives(1.0) - , diffusivity(PM_G2) + , diffusivity(cv::DIFF_PM_G2) , dthreshold(0.001f) , min_dthreshold(0.00001f) - , descriptor(cv::AKAZE::DESCRIPTOR_MLDB) + , descriptor(cv::DESCRIPTOR_MLDB) , descriptor_size(0) , descriptor_channels(3) , descriptor_pattern_size(10) @@ -67,12 +60,12 @@ struct AKAZEOptions { float soffset; ///< Base scale offset (sigma units) float derivative_factor; ///< Factor for the multiscale derivatives float sderivatives; ///< Smoothing factor for the derivatives - DIFFUSIVITY_TYPE diffusivity; ///< Diffusivity type + int diffusivity; ///< Diffusivity type float dthreshold; ///< Detector response threshold to accept point float min_dthreshold; ///< Minimum detector threshold to accept a point - cv::AKAZE::DESCRIPTOR_TYPE descriptor; ///< Type of descriptor + int descriptor; ///< Type of descriptor int descriptor_size; ///< Size of the descriptor in bits. 0->Full size int descriptor_channels; ///< Number of channels in the descriptor (1, 2, 3) int descriptor_pattern_size; ///< Actual patch size is 2*pattern_size*point.scale @@ -82,28 +75,4 @@ struct AKAZEOptions { int kcontrast_nbins; ///< Number of bins for the contrast factor histogram }; -/* ************************************************************************* */ -/// AKAZE nonlinear diffusion filtering evolution -struct TEvolution { - - TEvolution() { - etime = 0.0f; - esigma = 0.0f; - octave = 0; - sublevel = 0; - sigma_size = 0; - } - - 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 - size_t octave; // Image octave - size_t sublevel; // Image sublevel in each octave - size_t sigma_size; // Integer sigma. For computing the feature detector responses -}; \ No newline at end of file +#endif diff --git a/modules/features2d/src/kaze/AKAZEFeatures.cpp b/modules/features2d/src/kaze/AKAZEFeatures.cpp new file mode 100644 index 0000000..97222df --- /dev/null +++ b/modules/features2d/src/kaze/AKAZEFeatures.cpp @@ -0,0 +1,1880 @@ +/** + * @file AKAZEFeatures.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 "AKAZEFeatures.h" +#include "fed.h" +#include "nldiffusion_functions.h" +#include "utils.h" + +#include + +// Namespaces +using namespace std; +using namespace cv; +using namespace cv::details::kaze; + +/* ************************************************************************* */ +/** + * @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; + + if (options_.descriptor_size > 0 && options_.descriptor >= cv::DESCRIPTOR_MLDB_UPRIGHT) { + generateDescriptorSubsample(descriptorSamples_, descriptorBits_, options_.descriptor_size, + options_.descriptor_pattern_size, options_.descriptor_channels); + } + + Allocate_Memory_Evolution(); +} + +/* ************************************************************************* */ +/** + * @brief This method allocates the memory for the nonlinear diffusion evolution + */ +void AKAZEFeatures::Allocate_Memory_Evolution(void) { + + float rfactor = 0.0f; + 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.0f / 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; + } + + 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.Lsmooth = 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.5f*(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.0f; + ttime = evolution_[i].etime - evolution_[i - 1].etime; + naux = fed_tau_by_process_time(ttime, 1, 0.25f, 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 AKAZEFeatures::Create_Nonlinear_Scale_Space(const cv::Mat& img) +{ + CV_Assert(evolution_.size() > 0); + + // 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); + + // Allocate memory for the flow and step images + cv::Mat Lflow = cv::Mat::zeros(evolution_[0].Lt.rows, evolution_[0].Lt.cols, CV_32F); + cv::Mat Lstep = cv::Mat::zeros(evolution_[0].Lt.rows, evolution_[0].Lt.cols, CV_32F); + + // First compute the kcontrast factor + options_.kcontrast = compute_k_percentile(img, options_.kcontrast_percentile, 1.0f, options_.kcontrast_nbins, 0, 0); + + // 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.75f; + + // Allocate memory for the resized flow and step images + Lflow = cv::Mat::zeros(evolution_[i].Lt.rows, evolution_[i].Lt.cols, CV_32F); + Lstep = cv::Mat::zeros(evolution_[i].Lt.rows, evolution_[i].Lt.cols, CV_32F); + } + else { + evolution_[i - 1].Lt.copyTo(evolution_[i].Lt); + } + + gaussian_2D_convolution(evolution_[i].Lt, evolution_[i].Lsmooth, 0, 0, 1.0f); + + // 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 cv::DIFF_PM_G1: + pm_g1(evolution_[i].Lx, evolution_[i].Ly, Lflow, options_.kcontrast); + break; + case cv::DIFF_PM_G2: + pm_g2(evolution_[i].Lx, evolution_[i].Ly, Lflow, options_.kcontrast); + break; + case cv::DIFF_WEICKERT: + weickert_diffusivity(evolution_[i].Lx, evolution_[i].Ly, Lflow, options_.kcontrast); + break; + case cv::DIFF_CHARBONNIER: + charbonnier_diffusivity(evolution_[i].Lx, evolution_[i].Ly, Lflow, options_.kcontrast); + break; + default: + CV_Error(options_.diffusivity, "Diffusivity is not supported"); + break; + } + + // Perform FED n inner steps + for (int j = 0; j < nsteps_[i - 1]; j++) { + cv::details::kaze::nld_step_scalar(evolution_[i].Lt, Lflow, Lstep, tsteps_[i - 1][j]); + } + } + + 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) +{ + kpts.clear(); + Compute_Determinant_Hessian_Response(); + Find_Scale_Space_Extrema(kpts); + Do_Subpixel_Refinement(kpts); +} + +/* ************************************************************************* */ +class MultiscaleDerivativesAKAZEInvoker : public cv::ParallelLoopBody +{ +public: + explicit MultiscaleDerivativesAKAZEInvoker(std::vector& ev, const AKAZEOptions& opt) + : evolution_(&ev) + , options_(opt) + { + } + + void operator()(const cv::Range& range) const + { + std::vector& evolution = *evolution_; + + for (int i = range.start; i < range.end; i++) + { + 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_); + + 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_)); + } + } + +private: + std::vector* evolution_; + AKAZEOptions options_; +}; + +/* ************************************************************************* */ +/** + * @brief This method computes the multiscale derivatives for the nonlinear scale space + */ +void AKAZEFeatures::Compute_Multiscale_Derivatives(void) +{ + cv::parallel_for_(cv::Range(0, (int)evolution_.size()), + MultiscaleDerivativesAKAZEInvoker(evolution_, options_)); +} + +/* ************************************************************************* */ +/** + * @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(); + + for (size_t i = 0; i < evolution_.size(); i++) + { + 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) +{ + + 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 == cv::DESCRIPTOR_MLDB_UPRIGHT || options_.descriptor == cv::DESCRIPTOR_MLDB) { + smax = 10.0f*sqrtf(2.0f); + } + else if (options_.descriptor == cv::DESCRIPTOR_KAZE_UPRIGHT || options_.descriptor == cv::DESCRIPTOR_KAZE) { + smax = 12.0f*sqrtf(2.0f); + } + + 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 = (int)evolution_[i].octave; + point.class_id = (int)i; + ratio = pow(2.f, point.octave); + sigma_size_ = fRound(point.size / ratio); + point.pt.x = static_cast(jx); + point.pt.y = static_cast(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 = (int)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& pt = kpts_aux[i]; + for (size_t j = i + 1; j < kpts_aux.size(); j++) { + + // Compare response with the upper scale + if ((pt.class_id + 1) == kpts_aux[j].class_id) { + dist = sqrt(pow(pt.pt.x - kpts_aux[j].pt.x, 2) + pow(pt.pt.y - kpts_aux[j].pt.y, 2)); + if (dist <= pt.size) { + if (pt.response < kpts_aux[j].response) { + is_repeated = true; + break; + } + } + } + } + + if (is_repeated == false) + kpts.push_back(pt); + } +} + +/* ************************************************************************* */ +/** + * @brief This method performs subpixel refinement of the detected keypoints + * @param kpts Vector of detected keypoints + */ +void AKAZEFeatures::Do_Subpixel_Refinement(std::vector& kpts) +{ + 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); + + 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.5f)*(*(evolution_[kpts[i].class_id].Ldet.ptr(y)+x + 1) + - *(evolution_[kpts[i].class_id].Ldet.ptr(y)+x - 1)); + Dy = (0.5f)*(*(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.0f*(*(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.0f*(*(evolution_[kpts[i].class_id].Ldet.ptr(y)+x))); + + Dxy = (0.25f)*(*(evolution_[kpts[i].class_id].Ldet.ptr(y + 1) + x + 1) + + (*(evolution_[kpts[i].class_id].Ldet.ptr(y - 1) + x - 1))) + - (0.25f)*(*(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; + + cv::solve(A, b, dst, DECOMP_LU); + + if (fabs(*(dst.ptr(0))) <= 1.0f && fabs(*(dst.ptr(1))) <= 1.0f) { + kpts[i].pt.x = x + (*(dst.ptr(0))); + kpts[i].pt.y = y + (*(dst.ptr(1))); + kpts[i].pt.x *= powf(2.f, (float)evolution_[kpts[i].class_id].octave); + kpts[i].pt.y *= powf(2.f, (float)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.0f; + } + // Delete the point since its not stable + else { + kpts.erase(kpts.begin() + i); + i--; + } + } +} + +/* ************************************************************************* */ + +class SURF_Descriptor_Upright_64_Invoker : public cv::ParallelLoopBody +{ +public: + SURF_Descriptor_Upright_64_Invoker(std::vector& kpts, cv::Mat& desc, std::vector& evolution) + : keypoints_(&kpts) + , descriptors_(&desc) + , evolution_(&evolution) + { + } + + void operator() (const Range& range) const + { + for (int i = range.start; i < range.end; i++) + { + Get_SURF_Descriptor_Upright_64((*keypoints_)[i], descriptors_->ptr(i)); + } + } + + void Get_SURF_Descriptor_Upright_64(const cv::KeyPoint& kpt, float* desc) const; + +private: + std::vector* keypoints_; + cv::Mat* descriptors_; + std::vector* evolution_; +}; + +class SURF_Descriptor_64_Invoker : public cv::ParallelLoopBody +{ +public: + SURF_Descriptor_64_Invoker(std::vector& kpts, cv::Mat& desc, std::vector& evolution) + : keypoints_(&kpts) + , descriptors_(&desc) + , evolution_(&evolution) + { + } + + void operator()(const Range& range) const + { + for (int i = range.start; i < range.end; i++) + { + AKAZEFeatures::Compute_Main_Orientation((*keypoints_)[i], *evolution_); + Get_SURF_Descriptor_64((*keypoints_)[i], descriptors_->ptr(i)); + } + } + + void Get_SURF_Descriptor_64(const cv::KeyPoint& kpt, float* desc) const; + +private: + std::vector* keypoints_; + cv::Mat* descriptors_; + std::vector* evolution_; +}; + +class MSURF_Upright_Descriptor_64_Invoker : public cv::ParallelLoopBody +{ +public: + MSURF_Upright_Descriptor_64_Invoker(std::vector& kpts, cv::Mat& desc, std::vector& evolution) + : keypoints_(&kpts) + , descriptors_(&desc) + , evolution_(&evolution) + { + } + + void operator()(const Range& range) const + { + for (int i = range.start; i < range.end; i++) + { + Get_MSURF_Upright_Descriptor_64((*keypoints_)[i], descriptors_->ptr(i)); + } + } + + void Get_MSURF_Upright_Descriptor_64(const cv::KeyPoint& kpt, float* desc) const; + +private: + std::vector* keypoints_; + cv::Mat* descriptors_; + std::vector* evolution_; +}; + +class MSURF_Descriptor_64_Invoker : public cv::ParallelLoopBody +{ +public: + MSURF_Descriptor_64_Invoker(std::vector& kpts, cv::Mat& desc, std::vector& evolution) + : keypoints_(&kpts) + , descriptors_(&desc) + , evolution_(&evolution) + { + } + + void operator() (const Range& range) const + { + for (int i = range.start; i < range.end; i++) + { + AKAZEFeatures::Compute_Main_Orientation((*keypoints_)[i], *evolution_); + Get_MSURF_Descriptor_64((*keypoints_)[i], descriptors_->ptr(i)); + } + } + + void Get_MSURF_Descriptor_64(const cv::KeyPoint& kpt, float* desc) const; + +private: + std::vector* keypoints_; + cv::Mat* descriptors_; + std::vector* evolution_; +}; + +class Upright_MLDB_Full_Descriptor_Invoker : public cv::ParallelLoopBody +{ +public: + Upright_MLDB_Full_Descriptor_Invoker(std::vector& kpts, cv::Mat& desc, std::vector& evolution, AKAZEOptions& options) + : keypoints_(&kpts) + , descriptors_(&desc) + , evolution_(&evolution) + , options_(&options) + { + } + + void operator() (const Range& range) const + { + for (int i = range.start; i < range.end; i++) + { + Get_Upright_MLDB_Full_Descriptor((*keypoints_)[i], descriptors_->ptr(i)); + } + } + + void Get_Upright_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned char* desc) const; + +private: + std::vector* keypoints_; + cv::Mat* descriptors_; + std::vector* evolution_; + AKAZEOptions* options_; +}; + +class Upright_MLDB_Descriptor_Subset_Invoker : public cv::ParallelLoopBody +{ +public: + Upright_MLDB_Descriptor_Subset_Invoker(std::vector& kpts, + cv::Mat& desc, + std::vector& evolution, + AKAZEOptions& options, + cv::Mat descriptorSamples, + cv::Mat descriptorBits) + : keypoints_(&kpts) + , descriptors_(&desc) + , evolution_(&evolution) + , options_(&options) + , descriptorSamples_(descriptorSamples) + , descriptorBits_(descriptorBits) + { + } + + void operator() (const Range& range) const + { + for (int i = range.start; i < range.end; i++) + { + Get_Upright_MLDB_Descriptor_Subset((*keypoints_)[i], descriptors_->ptr(i)); + } + } + + void Get_Upright_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char* desc) const; + +private: + std::vector* keypoints_; + cv::Mat* descriptors_; + std::vector* evolution_; + AKAZEOptions* options_; + + cv::Mat descriptorSamples_; // List of positions in the grids to sample LDB bits from. + cv::Mat descriptorBits_; +}; + +class MLDB_Full_Descriptor_Invoker : public cv::ParallelLoopBody +{ +public: + MLDB_Full_Descriptor_Invoker(std::vector& kpts, cv::Mat& desc, std::vector& evolution, AKAZEOptions& options) + : keypoints_(&kpts) + , descriptors_(&desc) + , evolution_(&evolution) + , options_(&options) + { + } + + void operator() (const Range& range) const + { + for (int i = range.start; i < range.end; i++) + { + AKAZEFeatures::Compute_Main_Orientation((*keypoints_)[i], *evolution_); + Get_MLDB_Full_Descriptor((*keypoints_)[i], descriptors_->ptr(i)); + } + } + + void Get_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned char* desc) const; + +private: + std::vector* keypoints_; + cv::Mat* descriptors_; + std::vector* evolution_; + AKAZEOptions* options_; +}; + +class MLDB_Descriptor_Subset_Invoker : public cv::ParallelLoopBody +{ +public: + MLDB_Descriptor_Subset_Invoker(std::vector& kpts, + cv::Mat& desc, + std::vector& evolution, + AKAZEOptions& options, + cv::Mat descriptorSamples, + cv::Mat descriptorBits) + : keypoints_(&kpts) + , descriptors_(&desc) + , evolution_(&evolution) + , options_(&options) + , descriptorSamples_(descriptorSamples) + , descriptorBits_(descriptorBits) + { + } + + void operator() (const Range& range) const + { + for (int i = range.start; i < range.end; i++) + { + AKAZEFeatures::Compute_Main_Orientation((*keypoints_)[i], *evolution_); + Get_MLDB_Descriptor_Subset((*keypoints_)[i], descriptors_->ptr(i)); + } + } + + void Get_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char* desc) const; + +private: + std::vector* keypoints_; + cv::Mat* descriptors_; + std::vector* evolution_; + AKAZEOptions* options_; + + cv::Mat descriptorSamples_; // List of positions in the grids to sample LDB bits from. + cv::Mat descriptorBits_; +}; + +/** + * @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) +{ + for(size_t i = 0; i < kpts.size(); i++) + { + CV_Assert(0 <= kpts[i].class_id && kpts[i].class_id < static_cast(evolution_.size())); + } + + // Allocate memory for the matrix with the descriptors + if (options_.descriptor < cv::DESCRIPTOR_MLDB_UPRIGHT) { + desc = cv::Mat::zeros((int)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((int)kpts.size(), (int)ceil(t / 8.), CV_8UC1); + } + else { + // We use the random bit selection length binary descriptor + desc = cv::Mat::zeros((int)kpts.size(), (int)ceil(options_.descriptor_size / 8.), CV_8UC1); + } + } + + switch (options_.descriptor) + { + case cv::DESCRIPTOR_KAZE_UPRIGHT: // Upright descriptors, not invariant to rotation + { + cv::parallel_for_(cv::Range(0, (int)kpts.size()), MSURF_Upright_Descriptor_64_Invoker(kpts, desc, evolution_)); + } + break; + case cv::DESCRIPTOR_KAZE: + { + cv::parallel_for_(cv::Range(0, (int)kpts.size()), MSURF_Descriptor_64_Invoker(kpts, desc, evolution_)); + } + break; + case cv::DESCRIPTOR_MLDB_UPRIGHT: // Upright descriptors, not invariant to rotation + { + if (options_.descriptor_size == 0) + cv::parallel_for_(cv::Range(0, (int)kpts.size()), Upright_MLDB_Full_Descriptor_Invoker(kpts, desc, evolution_, options_)); + else + cv::parallel_for_(cv::Range(0, (int)kpts.size()), Upright_MLDB_Descriptor_Subset_Invoker(kpts, desc, evolution_, options_, descriptorSamples_, descriptorBits_)); + } + break; + case cv::DESCRIPTOR_MLDB: + { + if (options_.descriptor_size == 0) + cv::parallel_for_(cv::Range(0, (int)kpts.size()), MLDB_Full_Descriptor_Invoker(kpts, desc, evolution_, options_)); + else + cv::parallel_for_(cv::Range(0, (int)kpts.size()), MLDB_Descriptor_Subset_Invoker(kpts, desc, evolution_, options_, descriptorSamples_, descriptorBits_)); + } + break; + } +} + +/* ************************************************************************* */ +/** + * @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 AKAZEFeatures::Compute_Main_Orientation(cv::KeyPoint& kpt, const std::vector& evolution_) { + + 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 << evolution_[level].octave); + s = fRound(0.5f*kpt.size / ratio); + xf = kpt.pt.x / ratio; + yf = kpt.pt.y / ratio; + + // 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); + + 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] = getAngle(resX[idx], resY[idx]); + ++idx; + } + } + } + // Loop slides pi/3 window around feature point + for (ang1 = 0; ang1 < (float)(2.0 * CV_PI); ang1 += 0.15f) { + ang2 = (ang1 + (float)(CV_PI / 3.0) >(float)(2.0*CV_PI) ? ang1 - (float)(5.0*CV_PI / 3.0) : ang1 + (float)(CV_PI / 3.0)); + 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.0f*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 (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 MSURF_Upright_Descriptor_64_Invoker::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; + + // Subregion centers for the 4x4 gaussian weighting + float cx = -0.5f, cy = 0.5f; + + const std::vector& evolution = *evolution_; + + // 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 << kpt.octave); + scale = fRound(0.5f*kpt.size / ratio); + level = kpt.class_id; + yf = kpt.pt.y / ratio; + xf = kpt.pt.x / ratio; + + 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.0f; + cy = -0.5f; + + while (j < pattern_size) { + dx = dy = mdx = mdy = 0.0; + cy += 1.0f; + 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.50f*scale); + + y1 = (int)(sample_y - .5); + x1 = (int)(sample_x - .5); + + y2 = (int)(sample_y + .5); + x2 = (int)(sample_x + .5); + + 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.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + (1.0f - 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.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + (1.0f - 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 (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 MSURF_Descriptor_64_Invoker::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; + + // Subregion centers for the 4x4 gaussian weighting + float cx = -0.5f, cy = 0.5f; + + const std::vector& evolution = *evolution_; + + // 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 << kpt.octave); + scale = fRound(0.5f*kpt.size / ratio); + angle = kpt.angle; + level = kpt.class_id; + yf = kpt.pt.y / ratio; + xf = kpt.pt.x / ratio; + 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.0f; + cy = -0.5f; + + while (j < pattern_size) { + dx = dy = mdx = mdy = 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.5f*scale); + + y1 = fRound(sample_y - 0.5f); + x1 = fRound(sample_x - 0.5f); + + y2 = fRound(sample_y + 0.5f); + x2 = fRound(sample_x + 0.5f); + + 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.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + (1.0f - 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.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + (1.0f - 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 (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 Upright_MLDB_Full_Descriptor_Invoker::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; + + const AKAZEOptions & options = *options_; + const std::vector& evolution = *evolution_; + + // 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 << kpt.octave); + scale = fRound(0.5f*kpt.size / ratio); + level = kpt.class_id; + yf = kpt.pt.y / ratio; + xf = kpt.pt.x / ratio; + + // First 2x2 grid + pattern_size = options_->descriptor_pattern_size; + sample_step = pattern_size; + + 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_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 = static_cast(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++; + } + } + + //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_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++; + } + } + + //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++; + } + } +} + +/* ************************************************************************* */ +/** + * @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 MLDB_Full_Descriptor_Invoker::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; + + const AKAZEOptions & options = *options_; + const std::vector& evolution = *evolution_; + + // 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 << kpt.octave); + scale = fRound(0.5f*kpt.size / ratio); + angle = kpt.angle; + level = kpt.class_id; + yf = kpt.pt.y / ratio; + xf = kpt.pt.x / ratio; + co = cos(angle); + si = sin(angle); + + // First 2x2 grid + pattern_size = options.descriptor_pattern_size; + sample_step = pattern_size; + + 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 (float k = (float)i; k < i + sample_step; k++) { + for (float l = (float)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_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) { + 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 (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++; + } + } + } + + // Second 3x3 grid + sample_step = static_cast(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++; + } + } + + // 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 (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++; + } + } + } + + 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++; + } + } + } + + // 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++; + } + } + + // 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 (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++; + } + } + } +} + +/* ************************************************************************* */ +/** + * @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 MLDB_Descriptor_Subset_Invoker::Get_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char *desc) const { + + 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; + + const AKAZEOptions & options = *options_; + const std::vector& evolution = *evolution_; + + // Get the information from the keypoint + float ratio = (float)(1 << kpt.octave); + int scale = fRound(0.5f*kpt.size / ratio); + float angle = kpt.angle; + int level = kpt.class_id; + float yf = kpt.pt.y / ratio; + float xf = kpt.pt.x / ratio; + float co = cos(angle); + float si = sin(angle); + + // 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) = (int)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++) { + const 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*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 (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; + } + } + } + } + + *(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; + } + } + + // 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 Upright_MLDB_Descriptor_Subset_Invoker::Get_Upright_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char *desc) const { + + 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; + + const AKAZEOptions & options = *options_; + const std::vector& evolution = *evolution_; + + // Get the information from the keypoint + float ratio = (float)(1 << kpt.octave); + int scale = fRound(0.5f*kpt.size / ratio); + int level = kpt.class_id; + float yf = kpt.pt.y / ratio; + float xf = kpt.pt.x / ratio; + + // 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) = static_cast(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++) { + const 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 (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; + } + } + } + } + + *(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; + } + } + + // 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 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 (int i = 0, c = 0; i < 3; i++) { + int gdiv = i + 2; //grid divisions, per row + int gsz = gdiv*gdiv; + int psz = (int)ceil(2.f*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; + } + } + } + + srand(1024); + Mat_ comps = Mat_(nchannels * (int)ceil(nbits / (float)nchannels), 2); + comps = 1000; + + // Select some samples. A sample includes all channels + int count = 0; + int npicks = (int)ceil(nbits / (float)nchannels); + Mat_ samples(29, 3); + Mat_ fullcopy = fullM.clone(); + samples = -1; + + for (int i = 0; i < npicks; i++) { + int 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); + } + + sampleList = samples.rowRange(0, count).clone(); + comparisons = comps.rowRange(0, nbits).clone(); +} diff --git a/modules/features2d/src/kaze/AKAZEFeatures.h b/modules/features2d/src/kaze/AKAZEFeatures.h new file mode 100644 index 0000000..f8ce7a4 --- /dev/null +++ b/modules/features2d/src/kaze/AKAZEFeatures.h @@ -0,0 +1,62 @@ +/** + * @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 __OPENCV_FEATURES_2D_AKAZE_FEATURES_H__ +#define __OPENCV_FEATURES_2D_AKAZE_FEATURES_H__ + +/* ************************************************************************* */ +// Includes +#include "precomp.hpp" +#include "AKAZEConfig.h" +#include "TEvolution.h" + +/* ************************************************************************* */ +// AKAZE Class Declaration +class AKAZEFeatures { + +private: + + AKAZEOptions options_; ///< Configuration options for AKAZE + 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 + + /// 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_; + +public: + + /// Constructor with input arguments + AKAZEFeatures(const AKAZEOptions& options); + + /// Scale Space methods + void Allocate_Memory_Evolution(); + 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); + + /// Feature description methods + void Compute_Descriptors(std::vector& kpts, cv::Mat& desc); + static void Compute_Main_Orientation(cv::KeyPoint& kpt, const std::vector& evolution_); +}; + +/* ************************************************************************* */ +/// Inline functions +void generateDescriptorSubsample(cv::Mat& sampleList, cv::Mat& comparisons, + int nbits, int pattern_size, int nchannels); + +#endif diff --git a/modules/features2d/src/kaze/KAZEConfig.h b/modules/features2d/src/kaze/KAZEConfig.h index 988e247..21489a0 100644 --- a/modules/features2d/src/kaze/KAZEConfig.h +++ b/modules/features2d/src/kaze/KAZEConfig.h @@ -5,7 +5,8 @@ * @author Pablo F. Alcantarilla */ -#pragma once +#ifndef __OPENCV_FEATURES_2D_AKAZE_CONFIG_H__ +#define __OPENCV_FEATURES_2D_AKAZE_CONFIG_H__ // OpenCV Includes #include "precomp.hpp" @@ -15,14 +16,8 @@ struct KAZEOptions { - enum DIFFUSIVITY_TYPE { - PM_G1 = 0, - PM_G2 = 1, - WEICKERT = 2 - }; - KAZEOptions() - : diffusivity(PM_G2) + : diffusivity(cv::DIFF_PM_G2) , soffset(1.60f) , omax(4) @@ -33,20 +28,13 @@ struct KAZEOptions { , dthreshold(0.001f) , kcontrast(0.01f) , kcontrast_percentille(0.7f) - , kcontrast_bins(300) - - , use_fed(true) + , kcontrast_bins(300) , upright(false) , extended(false) - - , use_clipping_normalilzation(false) - , clipping_normalization_ratio(1.6f) - , clipping_normalization_niter(5) { } - DIFFUSIVITY_TYPE diffusivity; - + int diffusivity; float soffset; int omax; int nsublevels; @@ -57,27 +45,8 @@ struct KAZEOptions { float kcontrast; float kcontrast_percentille; int kcontrast_bins; - - bool use_fed; bool upright; bool extended; - - bool use_clipping_normalilzation; - float clipping_normalization_ratio; - int clipping_normalization_niter; }; -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/KAZEFeatures.cpp b/modules/features2d/src/kaze/KAZEFeatures.cpp index 634f68d..69e9a4b 100644 --- a/modules/features2d/src/kaze/KAZEFeatures.cpp +++ b/modules/features2d/src/kaze/KAZEFeatures.cpp @@ -22,22 +22,21 @@ */ #include "KAZEFeatures.h" +#include "utils.h" // Namespaces using namespace std; using namespace cv; using namespace cv::details::kaze; -//******************************************************************************* -//******************************************************************************* - +/* ************************************************************************* */ /** * @brief KAZE constructor with input options * @param options KAZE configuration options * @note The constructor allocates memory for the nonlinear scale space */ -KAZEFeatures::KAZEFeatures(KAZEOptions& _options) - : options(_options) +KAZEFeatures::KAZEFeatures(KAZEOptions& options) + : options_(options) { ncycles_ = 0; reordering_ = true; @@ -46,70 +45,48 @@ KAZEFeatures::KAZEFeatures(KAZEOptions& _options) Allocate_Memory_Evolution(); } -//******************************************************************************* -//******************************************************************************* - +/* ************************************************************************* */ /** * @brief This method allocates the memory for the nonlinear diffusion evolution */ void KAZEFeatures::Allocate_Memory_Evolution(void) { // Allocate the dimension of the matrices for the evolution - for (int i = 0; i <= options.omax - 1; i++) { - for (int j = 0; j <= options.nsublevels - 1; j++) { + for (int i = 0; i <= options_.omax - 1; i++) { + for (int j = 0; j <= options_.nsublevels - 1; j++) { TEvolution aux; - aux.Lx = cv::Mat::zeros(options.img_height, options.img_width, CV_32F); - aux.Ly = cv::Mat::zeros(options.img_height, options.img_width, CV_32F); - aux.Lxx = cv::Mat::zeros(options.img_height, options.img_width, CV_32F); - aux.Lxy = cv::Mat::zeros(options.img_height, options.img_width, CV_32F); - aux.Lyy = cv::Mat::zeros(options.img_height, options.img_width, CV_32F); - aux.Lflow = cv::Mat::zeros(options.img_height, options.img_width, CV_32F); - aux.Lt = cv::Mat::zeros(options.img_height, options.img_width, CV_32F); - aux.Lsmooth = cv::Mat::zeros(options.img_height, options.img_width, CV_32F); - aux.Lstep = cv::Mat::zeros(options.img_height, options.img_width, CV_32F); - aux.Ldet = cv::Mat::zeros(options.img_height, options.img_width, CV_32F); - aux.esigma = options.soffset*pow((float)2.0f, (float)(j) / (float)(options.nsublevels)+i); + aux.Lx = cv::Mat::zeros(options_.img_height, options_.img_width, CV_32F); + aux.Ly = cv::Mat::zeros(options_.img_height, options_.img_width, CV_32F); + aux.Lxx = cv::Mat::zeros(options_.img_height, options_.img_width, CV_32F); + aux.Lxy = cv::Mat::zeros(options_.img_height, options_.img_width, CV_32F); + aux.Lyy = cv::Mat::zeros(options_.img_height, options_.img_width, CV_32F); + aux.Lt = cv::Mat::zeros(options_.img_height, options_.img_width, CV_32F); + aux.Lsmooth = cv::Mat::zeros(options_.img_height, options_.img_width, CV_32F); + aux.Ldet = cv::Mat::zeros(options_.img_height, options_.img_width, CV_32F); + aux.esigma = options_.soffset*pow((float)2.0f, (float)(j) / (float)(options_.nsublevels)+i); aux.etime = 0.5f*(aux.esigma*aux.esigma); aux.sigma_size = fRound(aux.esigma); - aux.octave = (float)i; - aux.sublevel = (float)j; + aux.octave = i; + aux.sublevel = j; evolution_.push_back(aux); } } // Allocate memory for the FED number of cycles and time steps - if (options.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.25f, 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(options.img_width, options.img_height, CV_32F); // TODO? IS IT A BUG??? - Lty_ = Mat::zeros(options.img_height, options.img_width, CV_32F); - px_ = Mat::zeros(options.img_height, options.img_width, CV_32F); - py_ = Mat::zeros(options.img_height, options.img_width, CV_32F); - ax_ = Mat::zeros(options.img_height, options.img_width, CV_32F); - ay_ = Mat::zeros(options.img_height, options.img_width, CV_32F); - bx_ = Mat::zeros(options.img_height - 1, options.img_width, CV_32F); - by_ = Mat::zeros(options.img_height - 1, options.img_width, CV_32F); - qr_ = Mat::zeros(options.img_height - 1, options.img_width, CV_32F); - qc_ = Mat::zeros(options.img_height, options.img_width - 1, CV_32F); + 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.25f, 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 @@ -121,52 +98,47 @@ int KAZEFeatures::Create_Nonlinear_Scale_Space(const cv::Mat &img) // 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); - gaussian_2D_convolution(evolution_[0].Lt, evolution_[0].Lsmooth, 0, 0, options.sderivatives); + gaussian_2D_convolution(evolution_[0].Lt, evolution_[0].Lt, 0, 0, options_.soffset); + gaussian_2D_convolution(evolution_[0].Lt, evolution_[0].Lsmooth, 0, 0, options_.sderivatives); // Firstly compute the kcontrast factor - Compute_KContrast(evolution_[0].Lt, options.kcontrast_percentille); + Compute_KContrast(evolution_[0].Lt, options_.kcontrast_percentille); + + // Allocate memory for the flow and step images + cv::Mat Lflow = cv::Mat::zeros(evolution_[0].Lt.rows, evolution_[0].Lt.cols, CV_32F); + cv::Mat Lstep = cv::Mat::zeros(evolution_[0].Lt.rows, evolution_[0].Lt.cols, CV_32F); // 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, options.sderivatives); + gaussian_2D_convolution(evolution_[i - 1].Lt, evolution_[i].Lsmooth, 0, 0, options_.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 (options.diffusivity == KAZEOptions::PM_G1) { - pm_g1(evolution_[i].Lx, evolution_[i].Ly, evolution_[i].Lflow, options.kcontrast); + if (options_.diffusivity == cv::DIFF_PM_G1) { + pm_g1(evolution_[i].Lx, evolution_[i].Ly, Lflow, options_.kcontrast); } - else if (options.diffusivity == KAZEOptions::PM_G2) { - pm_g2(evolution_[i].Lx, evolution_[i].Ly, evolution_[i].Lflow, options.kcontrast); + else if (options_.diffusivity == cv::DIFF_PM_G2) { + pm_g2(evolution_[i].Lx, evolution_[i].Ly, Lflow, options_.kcontrast); } - else if (options.diffusivity == KAZEOptions::WEICKERT) { - weickert_diffusivity(evolution_[i].Lx, evolution_[i].Ly, evolution_[i].Lflow, options.kcontrast); + else if (options_.diffusivity == cv::DIFF_WEICKERT) { + weickert_diffusivity(evolution_[i].Lx, evolution_[i].Ly, Lflow, options_.kcontrast); } // Perform FED n inner steps - if (options.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); + for (int j = 0; j < nsteps_[i - 1]; j++) { + nld_step_scalar(evolution_[i].Lt, Lflow, Lstep, tsteps_[i - 1][j]); } } return 0; } -//************************************************************************************* -//************************************************************************************* - +/* ************************************************************************* */ /** * @brief This method computes the k contrast factor * @param img Input image @@ -174,38 +146,10 @@ int KAZEFeatures::Create_Nonlinear_Scale_Space(const cv::Mat &img) */ void KAZEFeatures::Compute_KContrast(const cv::Mat &img, const float &kpercentile) { - options.kcontrast = compute_k_percentile(img, kpercentile, options.sderivatives, options.kcontrast_bins, 0, 0); -} - -//************************************************************************************* -//************************************************************************************* - -/** - * @brief This method computes the multiscale derivatives for the nonlinear scale space - */ -void KAZEFeatures::Compute_Multiscale_Derivatives(void) -{ - // TODO: use cv::parallel_for_ - for (size_t i = 0; i < evolution_.size(); i++) - { - // 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)); - } + options_.kcontrast = compute_k_percentile(img, kpercentile, options_.sderivatives, options_.kcontrast_bins, 0, 0); } -//************************************************************************************* -//************************************************************************************* - +/* ************************************************************************* */ /** * @brief This method computes the feature detector response for the nonlinear scale space * @note We use the Hessian determinant as feature detector @@ -219,9 +163,9 @@ void KAZEFeatures::Compute_Detector_Response(void) for (size_t i = 0; i < evolution_.size(); i++) { - for (int ix = 0; ix < options.img_height; ix++) + for (int ix = 0; ix < options_.img_height; ix++) { - for (int jx = 0; jx < options.img_width; jx++) + for (int jx = 0; jx < options_.img_width; jx++) { lxx = *(evolution_[i].Lxx.ptr(ix)+jx); lxy = *(evolution_[i].Lxy.ptr(ix)+jx); @@ -232,9 +176,7 @@ void KAZEFeatures::Compute_Detector_Response(void) } } -//************************************************************************************* -//************************************************************************************* - +/* ************************************************************************* */ /** * @brief This method selects interesting keypoints through the nonlinear scale space * @param kpts Vector of keypoints @@ -242,27 +184,127 @@ void KAZEFeatures::Compute_Detector_Response(void) void KAZEFeatures::Feature_Detection(std::vector& kpts) { kpts.clear(); + Compute_Detector_Response(); + Determinant_Hessian(kpts); + Do_Subpixel_Refinement(kpts); +} - // Firstly compute the detector response for each pixel and scale level - Compute_Detector_Response(); +/* ************************************************************************* */ +class MultiscaleDerivativesKAZEInvoker : public cv::ParallelLoopBody +{ +public: + explicit MultiscaleDerivativesKAZEInvoker(std::vector& ev) : evolution_(&ev) + { + } + + void operator()(const cv::Range& range) const + { + std::vector& evolution = *evolution_; + for (int i = range.start; i < range.end; i++) + { + 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)); + } + } - // Find scale space extrema - Determinant_Hessian_Parallel(kpts); +private: + std::vector* evolution_; +}; - // Perform some subpixel refinement - Do_Subpixel_Refinement(kpts); +/* ************************************************************************* */ +/** + * @brief This method computes the multiscale derivatives for the nonlinear scale space + */ +void KAZEFeatures::Compute_Multiscale_Derivatives(void) +{ + cv::parallel_for_(cv::Range(0, (int)evolution_.size()), + MultiscaleDerivativesKAZEInvoker(evolution_)); } -//************************************************************************************* -//************************************************************************************* +/* ************************************************************************* */ +class FindExtremumKAZEInvoker : public cv::ParallelLoopBody +{ +public: + explicit FindExtremumKAZEInvoker(std::vector& ev, std::vector >& kpts_par, + const KAZEOptions& options) : evolution_(&ev), kpts_par_(&kpts_par), options_(options) + { + } + + void operator()(const cv::Range& range) const + { + std::vector& evolution = *evolution_; + std::vector >& kpts_par = *kpts_par_; + for (int i = range.start; i < range.end; i++) + { + float value = 0.0; + bool is_extremum = false; + + for (int ix = 1; ix < options_.img_height - 1; ix++) { + for (int jx = 1; jx < options_.img_width - 1; jx++) { + + is_extremum = false; + value = *(evolution[i].Ldet.ptr(ix)+jx); + + // Filter the points with the detector threshold + if (value > options_.dthreshold) { + if (value >= *(evolution[i].Ldet.ptr(ix)+jx - 1)) { + // First check on the same scale + if (check_maximum_neighbourhood(evolution[i].Ldet, 1, value, ix, jx, 1)) { + // Now check on the lower scale + if (check_maximum_neighbourhood(evolution[i - 1].Ldet, 1, value, ix, jx, 0)) { + // Now check on the upper scale + if (check_maximum_neighbourhood(evolution[i + 1].Ldet, 1, value, ix, jx, 0)) { + is_extremum = true; + } + } + } + } + } + + // Add the point of interest!! + if (is_extremum == true) { + cv::KeyPoint point; + point.pt.x = (float)jx; + point.pt.y = (float)ix; + point.response = fabs(value); + point.size = evolution[i].esigma; + point.octave = (int)evolution[i].octave; + point.class_id = i; + + // We use the angle field for the sublevel value + // Then, we will replace this angle field with the main orientation + point.angle = static_cast(evolution[i].sublevel); + kpts_par[i - 1].push_back(point); + } + } + } + } + } + +private: + std::vector* evolution_; + std::vector >* kpts_par_; + KAZEOptions options_; +}; + +/* ************************************************************************* */ /** * @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 KAZEFeatures::Determinant_Hessian_Parallel(std::vector& kpts) +void KAZEFeatures::Determinant_Hessian(std::vector& kpts) { int level = 0; float dist = 0.0, smax = 3.0; @@ -283,10 +325,8 @@ void KAZEFeatures::Determinant_Hessian_Parallel(std::vector& kpts) kpts_par_.push_back(aux); } - // TODO: Use cv::parallel_for_ - for (int i = 1; i < (int)evolution_.size() - 1; i++) { - Find_Extremum_Threading(i); - } + cv::parallel_for_(cv::Range(1, (int)evolution_.size()-1), + FindExtremumKAZEInvoker(evolution_, kpts_par_, options_)); // Now fill the vector of keypoints!!! for (int i = 0; i < (int)kpts_par_.size(); i++) { @@ -343,63 +383,7 @@ void KAZEFeatures::Determinant_Hessian_Parallel(std::vector& kpts) } } -//************************************************************************************* -//************************************************************************************* - -/** - * @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 KAZEFeatures::Find_Extremum_Threading(const int& level) { - - float value = 0.0; - bool is_extremum = false; - - for (int ix = 1; ix < options.img_height - 1; ix++) { - for (int jx = 1; jx < options.img_width - 1; jx++) { - - is_extremum = false; - value = *(evolution_[level].Ldet.ptr(ix)+jx); - - // Filter the points with the detector threshold - if (value > options.dthreshold) { - 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 = (float)jx; - point.pt.y = (float)ix; - point.response = fabs(value); - point.size = evolution_[level].esigma; - point.octave = (int)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 @@ -475,10 +459,10 @@ void KAZEFeatures::Do_Subpixel_Refinement(std::vector &kpts) { if (fabs(*(dst.ptr(0))) <= 1.0f && fabs(*(dst.ptr(1))) <= 1.0f && fabs(*(dst.ptr(2))) <= 1.0f) { 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)(options.nsublevels)); + dsc = kpts_[i].octave + (kpts_[i].angle + *(dst.ptr(2))) / ((float)(options_.nsublevels)); // In OpenCV the size of a keypoint is the diameter!! - kpts_[i].size = 2.0f*options.soffset*pow((float)2.0f, dsc); + kpts_[i].size = 2.0f*options_.soffset*pow((float)2.0f, dsc); kpts_[i].angle = 0.0; } // Set the points to be deleted after the for loop @@ -497,17 +481,15 @@ void KAZEFeatures::Do_Subpixel_Refinement(std::vector &kpts) { } } -//************************************************************************************* -//************************************************************************************* - +/* ************************************************************************* */ class KAZE_Descriptor_Invoker : public cv::ParallelLoopBody { public: - KAZE_Descriptor_Invoker(std::vector &kpts, cv::Mat &desc, std::vector& evolution, const KAZEOptions& _options) - : _kpts(&kpts) - , _desc(&desc) - , _evolution(&evolution) - , options(_options) + KAZE_Descriptor_Invoker(std::vector &kpts, cv::Mat &desc, std::vector& evolution, const KAZEOptions& options) + : kpts_(&kpts) + , desc_(&desc) + , evolution_(&evolution) + , options_(options) { } @@ -517,26 +499,26 @@ public: void operator() (const cv::Range& range) const { - std::vector &kpts = *_kpts; - cv::Mat &desc = *_desc; - std::vector &evolution = *_evolution; + std::vector &kpts = *kpts_; + cv::Mat &desc = *desc_; + std::vector &evolution = *evolution_; for (int i = range.start; i < range.end; i++) { kpts[i].angle = 0.0; - if (options.upright) + if (options_.upright) { kpts[i].angle = 0.0; - if (options.extended) + if (options_.extended) Get_KAZE_Upright_Descriptor_128(kpts[i], desc.ptr((int)i)); else Get_KAZE_Upright_Descriptor_64(kpts[i], desc.ptr((int)i)); } else { - KAZEFeatures::Compute_Main_Orientation(kpts[i], evolution, options); + KAZEFeatures::Compute_Main_Orientation(kpts[i], evolution, options_); - if (options.extended) + if (options_.extended) Get_KAZE_Descriptor_128(kpts[i], desc.ptr((int)i)); else Get_KAZE_Descriptor_64(kpts[i], desc.ptr((int)i)); @@ -549,12 +531,13 @@ private: void Get_KAZE_Upright_Descriptor_128(const cv::KeyPoint& kpt, float* desc) const; void Get_KAZE_Descriptor_128(const cv::KeyPoint& kpt, float *desc) const; - std::vector * _kpts; - cv::Mat * _desc; - std::vector * _evolution; - KAZEOptions options; + std::vector * kpts_; + cv::Mat * desc_; + std::vector * evolution_; + KAZEOptions options_; }; +/* ************************************************************************* */ /** * @brief This method computes the set of descriptors through the nonlinear scale space * @param kpts Vector of keypoints @@ -562,20 +545,23 @@ private: */ void KAZEFeatures::Feature_Description(std::vector &kpts, cv::Mat &desc) { + for(size_t i = 0; i < kpts.size(); i++) + { + CV_Assert(0 <= kpts[i].class_id && kpts[i].class_id < static_cast(evolution_.size())); + } + // Allocate memory for the matrix of descriptors - if (options.extended == true) { + if (options_.extended == true) { desc = Mat::zeros((int)kpts.size(), 128, CV_32FC1); } else { desc = Mat::zeros((int)kpts.size(), 64, CV_32FC1); } - cv::parallel_for_(cv::Range(0, (int)kpts.size()), KAZE_Descriptor_Invoker(kpts, desc, evolution_, options)); + cv::parallel_for_(cv::Range(0, (int)kpts.size()), KAZE_Descriptor_Invoker(kpts, desc, evolution_, options_)); } -//************************************************************************************* -//************************************************************************************* - +/* ************************************************************************* */ /** * @brief This method computes the main orientation for a given keypoint * @param kpt Input keypoint @@ -651,9 +637,7 @@ void KAZEFeatures::Compute_Main_Orientation(cv::KeyPoint &kpt, const std::vector } } -//************************************************************************************* -//************************************************************************************* - +/* ************************************************************************* */ /** * @brief This method computes the upright descriptor (not rotation invariant) of * the provided keypoint @@ -673,7 +657,7 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Upright_Descriptor_64(const cv::KeyPoint 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; - std::vector& evolution_ = *_evolution; + std::vector& evolution = *evolution_; // Subregion centers for the 4x4 gaussian weighting float cx = -0.5f, cy = 0.5f; @@ -724,26 +708,26 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Upright_Descriptor_64(const cv::KeyPoint y1 = (int)(sample_y - 0.5f); x1 = (int)(sample_x - 0.5f); - checkDescriptorLimits(x1, y1, options.img_width, options.img_height); + checkDescriptorLimits(x1, y1, options_.img_width, options_.img_height); y2 = (int)(sample_y + 0.5f); x2 = (int)(sample_x + 0.5f); - checkDescriptorLimits(x2, y2, options.img_width, options.img_height); + checkDescriptorLimits(x2, y2, options_.img_width, options_.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); + 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.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + (1.0f - 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); + 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.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + (1.0f - fx)*fy*res3 + fx*fy*res4; rx = gauss_s1*rx; @@ -779,15 +763,9 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Upright_Descriptor_64(const cv::KeyPoint for (i = 0; i < dsize; i++) { desc[i] /= len; } - - if (options.use_clipping_normalilzation) { - clippingDescriptor(desc, dsize, options.clipping_normalization_niter, options.clipping_normalization_ratio); - } } -//************************************************************************************* -//************************************************************************************* - +/* ************************************************************************* */ /** * @brief This method computes the descriptor of the provided keypoint given the * main orientation of the keypoint @@ -807,7 +785,7 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Descriptor_64(const cv::KeyPoint &kpt, fl int kx = 0, ky = 0, i = 0, j = 0, dcount = 0; int dsize = 0, scale = 0, level = 0; - std::vector& evolution_ = *_evolution; + std::vector& evolution = *evolution_; // Subregion centers for the 4x4 gaussian weighting float cx = -0.5f, cy = 0.5f; @@ -862,26 +840,26 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Descriptor_64(const cv::KeyPoint &kpt, fl y1 = fRound(sample_y - 0.5f); x1 = fRound(sample_x - 0.5f); - checkDescriptorLimits(x1, y1, options.img_width, options.img_height); + checkDescriptorLimits(x1, y1, options_.img_width, options_.img_height); y2 = (int)(sample_y + 0.5f); x2 = (int)(sample_x + 0.5f); - checkDescriptorLimits(x2, y2, options.img_width, options.img_height); + checkDescriptorLimits(x2, y2, options_.img_width, options_.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); + 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.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + (1.0f - 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); + 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.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + (1.0f - fx)*fy*res3 + fx*fy*res4; // Get the x and y derivatives on the rotated axis @@ -914,15 +892,9 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Descriptor_64(const cv::KeyPoint &kpt, fl for (i = 0; i < dsize; i++) { desc[i] /= len; } - - if (options.use_clipping_normalilzation) { - clippingDescriptor(desc, dsize, options.clipping_normalization_niter, options.clipping_normalization_ratio); - } } -//************************************************************************************* -//************************************************************************************* - +/* ************************************************************************* */ /** * @brief This method computes the extended upright descriptor (not rotation invariant) of * the provided keypoint @@ -947,7 +919,7 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Upright_Descriptor_128(const cv::KeyPoint // Subregion centers for the 4x4 gaussian weighting float cx = -0.5f, cy = 0.5f; - std::vector& evolution_ = *_evolution; + std::vector& evolution = *evolution_; // Set the descriptor size and the sample and pattern sizes dsize = 128; @@ -998,26 +970,26 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Upright_Descriptor_128(const cv::KeyPoint y1 = (int)(sample_y - 0.5f); x1 = (int)(sample_x - 0.5f); - checkDescriptorLimits(x1, y1, options.img_width, options.img_height); + checkDescriptorLimits(x1, y1, options_.img_width, options_.img_height); y2 = (int)(sample_y + 0.5f); x2 = (int)(sample_x + 0.5f); - checkDescriptorLimits(x2, y2, options.img_width, options.img_height); + checkDescriptorLimits(x2, y2, options_.img_width, options_.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); + 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.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + (1.0f - 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); + 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.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + (1.0f - fx)*fy*res3 + fx*fy*res4; rx = gauss_s1*rx; @@ -1072,15 +1044,9 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Upright_Descriptor_128(const cv::KeyPoint for (i = 0; i < dsize; i++) { desc[i] /= len; } - - if (options.use_clipping_normalilzation) { - clippingDescriptor(desc, dsize, options.clipping_normalization_niter, options.clipping_normalization_ratio); - } } -//************************************************************************************* -//************************************************************************************* - +/* ************************************************************************* */ /** * @brief This method computes the extended G-SURF descriptor of the provided keypoint * given the main orientation of the keypoint @@ -1102,7 +1068,7 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Descriptor_128(const cv::KeyPoint &kpt, f int kx = 0, ky = 0, i = 0, j = 0, dcount = 0; int dsize = 0, scale = 0, level = 0; - std::vector& evolution_ = *_evolution; + std::vector& evolution = *evolution_; // Subregion centers for the 4x4 gaussian weighting float cx = -0.5f, cy = 0.5f; @@ -1160,26 +1126,26 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Descriptor_128(const cv::KeyPoint &kpt, f y1 = fRound(sample_y - 0.5f); x1 = fRound(sample_x - 0.5f); - checkDescriptorLimits(x1, y1, options.img_width, options.img_height); + checkDescriptorLimits(x1, y1, options_.img_width, options_.img_height); y2 = (int)(sample_y + 0.5f); x2 = (int)(sample_x + 0.5f); - checkDescriptorLimits(x2, y2, options.img_width, options.img_height); + checkDescriptorLimits(x2, y2, options_.img_width, options_.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); + 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.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + (1.0f - 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); + 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.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + (1.0f - fx)*fy*res3 + fx*fy*res4; // Get the x and y derivatives on the rotated axis @@ -1235,298 +1201,4 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Descriptor_128(const cv::KeyPoint &kpt, f for (i = 0; i < dsize; i++) { desc[i] /= len; } - - if (options.use_clipping_normalilzation) { - clippingDescriptor(desc, dsize, options.clipping_normalization_niter, options.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 arbitrarily large - */ -void KAZEFeatures::AOS_Step_Scalar(cv::Mat &Ld, const cv::Mat &Ldprev, const cv::Mat &c, const float& stepsize) { - - AOS_Rows(Ldprev, c, stepsize); - AOS_Columns(Ldprev, c, stepsize); - - Ld = 0.5f*(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 KAZEFeatures::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.0f + 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 KAZEFeatures::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.0f + 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 KAZEFeatures::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 (float)CV_PI - atan(-y / x); - } - - if (x < 0 && y < 0) { - return (float)CV_PI + atan(y / x); - } - - if (x >= 0 && y < 0) { - return 2.0f * (float)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 / sqrtf(static_cast(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/KAZEFeatures.h b/modules/features2d/src/kaze/KAZEFeatures.h index 81509c4..b62f948 100644 --- a/modules/features2d/src/kaze/KAZEFeatures.h +++ b/modules/features2d/src/kaze/KAZEFeatures.h @@ -7,84 +7,53 @@ * @author Pablo F. Alcantarilla */ -#ifndef KAZE_H_ -#define KAZE_H_ - -//************************************************************************************* -//************************************************************************************* +#ifndef __OPENCV_FEATURES_2D_KAZE_FEATURES_H__ +#define __OPENCV_FEATURES_2D_KAZE_FEATURES_H__ +/* ************************************************************************* */ // Includes #include "KAZEConfig.h" #include "nldiffusion_functions.h" #include "fed.h" +#include "TEvolution.h" -//************************************************************************************* -//************************************************************************************* - +/* ************************************************************************* */ // KAZE Class Declaration class KAZEFeatures { private: - KAZEOptions options; - - // Parameters of the Nonlinear diffusion class - std::vector evolution_; // Vector of nonlinear diffusion evolution + /// Parameters of the Nonlinear diffusion class + KAZEOptions options_; ///< Configuration options for KAZE + std::vector evolution_; ///< Vector of nonlinear diffusion evolution - // Vector of keypoint vectors for finding extrema in multiple threads + /// 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 - - // Some auxiliary variables used in the AOS step - cv::Mat Ltx_, Lty_, px_, py_, ax_, ay_, bx_, by_, qr_, qc_; + /// 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 public: - // Constructor + /// Constructor KAZEFeatures(KAZEOptions& options); - // Public methods for KAZE interface + /// 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); - static void Compute_Main_Orientation(cv::KeyPoint& kpt, const std::vector& evolution_, const KAZEOptions& options); -private: - - // Feature Detection Methods + /// 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 Determinant_Hessian(std::vector& kpts); void Do_Subpixel_Refinement(std::vector& kpts); - - // 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); - }; -//************************************************************************************* -//************************************************************************************* - -// 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_ +#endif diff --git a/modules/features2d/src/kaze/TEvolution.h b/modules/features2d/src/kaze/TEvolution.h new file mode 100644 index 0000000..bc7654e --- /dev/null +++ b/modules/features2d/src/kaze/TEvolution.h @@ -0,0 +1,35 @@ +/** + * @file TEvolution.h + * @brief Header file with the declaration of the TEvolution struct + * @date Jun 02, 2014 + * @author Pablo F. Alcantarilla + */ + +#ifndef __OPENCV_FEATURES_2D_TEVOLUTION_H__ +#define __OPENCV_FEATURES_2D_TEVOLUTION_H__ + +/* ************************************************************************* */ +/// KAZE/A-KAZE nonlinear diffusion filtering evolution +struct TEvolution { + + TEvolution() { + etime = 0.0f; + esigma = 0.0f; + octave = 0; + sublevel = 0; + sigma_size = 0; + } + + cv::Mat Lx, Ly; ///< First order spatial derivatives + cv::Mat Lxx, Lxy, Lyy; ///< Second order spatial derivatives + cv::Mat Lt; ///< Evolution image + cv::Mat Lsmooth; ///< Smoothed image + 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 esigma. For computing the feature detector responses +}; + +#endif diff --git a/modules/features2d/src/kaze/fed.h b/modules/features2d/src/kaze/fed.h index c313b81..0cfa400 100644 --- a/modules/features2d/src/kaze/fed.h +++ b/modules/features2d/src/kaze/fed.h @@ -1,5 +1,5 @@ -#ifndef FED_H -#define FED_H +#ifndef __OPENCV_FEATURES_2D_FED_H__ +#define __OPENCV_FEATURES_2D_FED_H__ //****************************************************************************** //****************************************************************************** @@ -22,4 +22,4 @@ bool fed_is_prime_internal(const int& number); //************************************************************************************* //************************************************************************************* -#endif // FED_H +#endif // __OPENCV_FEATURES_2D_FED_H__ diff --git a/modules/features2d/src/kaze/nldiffusion_functions.h b/modules/features2d/src/kaze/nldiffusion_functions.h index 773f7e4..5c161a6 100644 --- a/modules/features2d/src/kaze/nldiffusion_functions.h +++ b/modules/features2d/src/kaze/nldiffusion_functions.h @@ -8,8 +8,8 @@ * @author Pablo F. Alcantarilla */ -#ifndef KAZE_NLDIFFUSION_FUNCTIONS_H -#define KAZE_NLDIFFUSION_FUNCTIONS_H +#ifndef __OPENCV_FEATURES_2D_NLDIFFUSION_FUNCTIONS_H__ +#define __OPENCV_FEATURES_2D_NLDIFFUSION_FUNCTIONS_H__ /* ************************************************************************* */ // Includes diff --git a/modules/features2d/src/kaze/utils.h b/modules/features2d/src/kaze/utils.h new file mode 100644 index 0000000..13ae352 --- /dev/null +++ b/modules/features2d/src/kaze/utils.h @@ -0,0 +1,77 @@ +#ifndef __OPENCV_FEATURES_2D_KAZE_UTILS_H__ +#define __OPENCV_FEATURES_2D_KAZE_UTILS_H__ + +/* ************************************************************************* */ +/** + * @brief This function computes the angle from the vector given by (X Y). From 0 to 2*Pi + */ +inline float getAngle(float x, float y) { + + if (x >= 0 && y >= 0) { + return atanf(y / x); + } + + if (x < 0 && y >= 0) { + return static_cast(CV_PI)-atanf(-y / x); + } + + if (x < 0 && y < 0) { + return static_cast(CV_PI)+atanf(y / x); + } + + if (x >= 0 && y < 0) { + return static_cast(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 checkDescriptorLimits(int &x, int &y, int width, 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); +} + +#endif diff --git a/modules/features2d/test/test_descriptors_regression.cpp b/modules/features2d/test/test_descriptors_regression.cpp index 281df24..3c6ae97 100644 --- a/modules/features2d/test/test_descriptors_regression.cpp +++ b/modules/features2d/test/test_descriptors_regression.cpp @@ -101,8 +101,14 @@ public: typedef typename Distance::ResultType DistanceType; CV_DescriptorExtractorTest( const string _name, DistanceType _maxDist, const Ptr& _dextractor, - Distance d = Distance() ): - name(_name), maxDist(_maxDist), dextractor(_dextractor), distance(d) {} + Distance d = Distance(), Ptr _detector = Ptr()): + name(_name), maxDist(_maxDist), dextractor(_dextractor), distance(d) , detector(_detector) {} + + ~CV_DescriptorExtractorTest() + { + if(!detector.empty()) + detector.release(); + } protected: virtual void createDescriptorExtractor() {} @@ -189,7 +195,6 @@ protected: // Read the test image. string imgFilename = string(ts->get_data_path()) + FEATURES2D_DIR + "/" + IMAGE_FILENAME; - Mat img = imread( imgFilename ); if( img.empty() ) { @@ -197,13 +202,15 @@ protected: ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_TEST_DATA ); return; } - vector keypoints; FileStorage fs( string(ts->get_data_path()) + FEATURES2D_DIR + "/keypoints.xml.gz", FileStorage::READ ); - if( fs.isOpened() ) - { + if(!detector.empty()) { + detector->detect(img, keypoints); + } else { read( fs.getFirstTopLevelNode(), keypoints ); - + } + if(!keypoints.empty()) + { Mat calcDescriptors; double t = (double)getTickCount(); dextractor->compute( img, keypoints, calcDescriptors ); @@ -244,7 +251,7 @@ protected: } } } - else + if(!fs.isOpened()) { ts->printf( cvtest::TS::LOG, "Compute and write keypoints.\n" ); fs.open( string(ts->get_data_path()) + FEATURES2D_DIR + "/keypoints.xml.gz", FileStorage::WRITE ); @@ -295,6 +302,7 @@ protected: const DistanceType maxDist; Ptr dextractor; Distance distance; + Ptr detector; private: CV_DescriptorExtractorTest& operator=(const CV_DescriptorExtractorTest&) { return *this; } @@ -340,3 +348,19 @@ TEST( Features2d_DescriptorExtractor_OpponentBRIEF, regression ) DescriptorExtractor::create("OpponentBRIEF") ); test.safe_run(); } + +TEST( Features2d_DescriptorExtractor_KAZE, regression ) +{ + CV_DescriptorExtractorTest< L2 > test( "descriptor-kaze", 0.03f, + DescriptorExtractor::create("KAZE"), + L2(), FeatureDetector::create("KAZE")); + test.safe_run(); +} + +TEST( Features2d_DescriptorExtractor_AKAZE, regression ) +{ + CV_DescriptorExtractorTest test( "descriptor-akaze", (CV_DescriptorExtractorTest::DistanceType)12.f, + DescriptorExtractor::create("AKAZE"), + Hamming(), FeatureDetector::create("AKAZE")); + test.safe_run(); +} diff --git a/modules/features2d/test/test_detectors_regression.cpp b/modules/features2d/test/test_detectors_regression.cpp index 8f34913..25e2c7f 100644 --- a/modules/features2d/test/test_detectors_regression.cpp +++ b/modules/features2d/test/test_detectors_regression.cpp @@ -289,6 +289,18 @@ TEST( Features2d_Detector_ORB, regression ) test.safe_run(); } +TEST( Features2d_Detector_KAZE, regression ) +{ + CV_FeatureDetectorTest test( "detector-kaze", FeatureDetector::create("KAZE") ); + test.safe_run(); +} + +TEST( Features2d_Detector_AKAZE, regression ) +{ + CV_FeatureDetectorTest test( "detector-akaze", FeatureDetector::create("AKAZE") ); + test.safe_run(); +} + TEST( Features2d_Detector_GridFAST, regression ) { CV_FeatureDetectorTest test( "detector-grid-fast", FeatureDetector::create("GridFAST") ); diff --git a/modules/features2d/test/test_keypoints.cpp b/modules/features2d/test/test_keypoints.cpp index 6f9a7e1..a9f30b1 100644 --- a/modules/features2d/test/test_keypoints.cpp +++ b/modules/features2d/test/test_keypoints.cpp @@ -167,19 +167,17 @@ TEST(Features2d_Detector_Keypoints_Dense, validation) test.safe_run(); } -// FIXIT #2807 Crash on Windows 7 x64 MSVS 2012, Linux Fedora 19 x64 with GCC 4.8.2, Linux Ubuntu 14.04 LTS x64 with GCC 4.8.2 -TEST(Features2d_Detector_Keypoints_KAZE, DISABLED_validation) +TEST(Features2d_Detector_Keypoints_KAZE, validation) { CV_FeatureDetectorKeypointsTest test(Algorithm::create("Feature2D.KAZE")); test.safe_run(); } -// FIXIT #2807 Crash on Windows 7 x64 MSVS 2012, Linux Fedora 19 x64 with GCC 4.8.2, Linux Ubuntu 14.04 LTS x64 with GCC 4.8.2 -TEST(Features2d_Detector_Keypoints_AKAZE, DISABLED_validation) +TEST(Features2d_Detector_Keypoints_AKAZE, validation) { - CV_FeatureDetectorKeypointsTest test_kaze(cv::Ptr(new cv::AKAZE(cv::AKAZE::DESCRIPTOR_KAZE))); + CV_FeatureDetectorKeypointsTest test_kaze(cv::Ptr(new cv::AKAZE(cv::DESCRIPTOR_KAZE))); test_kaze.safe_run(); - CV_FeatureDetectorKeypointsTest test_mldb(cv::Ptr(new cv::AKAZE(cv::AKAZE::DESCRIPTOR_MLDB))); + CV_FeatureDetectorKeypointsTest test_mldb(cv::Ptr(new cv::AKAZE(cv::DESCRIPTOR_MLDB))); test_mldb.safe_run(); } diff --git a/modules/features2d/test/test_rotation_and_scale_invariance.cpp b/modules/features2d/test/test_rotation_and_scale_invariance.cpp index 69ba6f0..fa6a136 100644 --- a/modules/features2d/test/test_rotation_and_scale_invariance.cpp +++ b/modules/features2d/test/test_rotation_and_scale_invariance.cpp @@ -652,8 +652,7 @@ TEST(Features2d_ScaleInvariance_Detector_BRISK, regression) test.safe_run(); } -// FIXIT #2807 Crash on Windows 7 x64 MSVS 2012, Linux Fedora 19 x64 with GCC 4.8.2, Linux Ubuntu 14.04 LTS x64 with GCC 4.8.2 -TEST(Features2d_ScaleInvariance_Detector_KAZE, DISABLED_regression) +TEST(Features2d_ScaleInvariance_Detector_KAZE, regression) { DetectorScaleInvarianceTest test(Algorithm::create("Feature2D.KAZE"), 0.08f, @@ -661,8 +660,7 @@ TEST(Features2d_ScaleInvariance_Detector_KAZE, DISABLED_regression) test.safe_run(); } -// FIXIT #2807 Crash on Windows 7 x64 MSVS 2012, Linux Fedora 19 x64 with GCC 4.8.2, Linux Ubuntu 14.04 LTS x64 with GCC 4.8.2 -TEST(Features2d_ScaleInvariance_Detector_AKAZE, DISABLED_regression) +TEST(Features2d_ScaleInvariance_Detector_AKAZE, regression) { DetectorScaleInvarianceTest test(Algorithm::create("Feature2D.AKAZE"), 0.08f, diff --git a/samples/cpp/H1to3p.xml b/samples/cpp/H1to3p.xml new file mode 100755 index 0000000..d8b7087 --- /dev/null +++ b/samples/cpp/H1to3p.xml @@ -0,0 +1,11 @@ + + + + 3 + 3 +
d
+ + 7.6285898e-01 -2.9922929e-01 2.2567123e+02 + 3.3443473e-01 1.0143901e+00 -7.6999973e+01 + 3.4663091e-04 -1.4364524e-05 1.0000000e+00 + diff --git a/samples/cpp/graf1.png b/samples/cpp/graf1.png new file mode 100755 index 0000000..67e3473 Binary files /dev/null and b/samples/cpp/graf1.png differ diff --git a/samples/cpp/graf3.png b/samples/cpp/graf3.png new file mode 100755 index 0000000..3bcc609 Binary files /dev/null and b/samples/cpp/graf3.png differ diff --git a/samples/cpp/tutorial_code/features2D/AKAZE_match.cpp b/samples/cpp/tutorial_code/features2D/AKAZE_match.cpp new file mode 100755 index 0000000..894627e --- /dev/null +++ b/samples/cpp/tutorial_code/features2D/AKAZE_match.cpp @@ -0,0 +1,79 @@ +#include +#include +#include +#include +#include + +using namespace std; +using namespace cv; + +const float inlier_threshold = 2.5f; // Distance threshold to identify inliers +const float nn_match_ratio = 0.8f; // Nearest neighbor matching ratio + +int main(void) +{ + Mat img1 = imread("graf1.png", IMREAD_GRAYSCALE); + Mat img2 = imread("graf3.png", IMREAD_GRAYSCALE); + + Mat homography; + FileStorage fs("H1to3p.xml", FileStorage::READ); + fs.getFirstTopLevelNode() >> homography; + + vector kpts1, kpts2; + Mat desc1, desc2; + + AKAZE akaze; + akaze(img1, noArray(), kpts1, desc1); + akaze(img2, noArray(), kpts2, desc2); + + BFMatcher matcher(NORM_HAMMING); + vector< vector > nn_matches; + matcher.knnMatch(desc1, desc2, nn_matches, 2); + + vector matched1, matched2, inliers1, inliers2; + vector good_matches; + for(size_t i = 0; i < nn_matches.size(); i++) { + DMatch first = nn_matches[i][0]; + float dist1 = nn_matches[i][0].distance; + float dist2 = nn_matches[i][1].distance; + + if(dist1 < nn_match_ratio * dist2) { + matched1.push_back(kpts1[first.queryIdx]); + matched2.push_back(kpts2[first.trainIdx]); + } + } + + for(unsigned i = 0; i < matched1.size(); i++) { + Mat col = Mat::ones(3, 1, CV_64F); + col.at(0) = matched1[i].pt.x; + col.at(1) = matched1[i].pt.y; + + col = homography * col; + col /= col.at(2); + double dist = sqrt( pow(col.at(0) - matched2[i].pt.x, 2) + + pow(col.at(1) - matched2[i].pt.y, 2)); + + if(dist < inlier_threshold) { + int new_i = static_cast(inliers1.size()); + inliers1.push_back(matched1[i]); + inliers2.push_back(matched2[i]); + good_matches.push_back(DMatch(new_i, new_i, 0)); + } + } + + Mat res; + drawMatches(img1, inliers1, img2, inliers2, good_matches, res); + imwrite("res.png", res); + + double inlier_ratio = inliers1.size() * 1.0 / matched1.size(); + cout << "A-KAZE Matching Results" << endl; + cout << "*******************************" << endl; + cout << "# Keypoints 1: \t" << kpts1.size() << endl; + cout << "# Keypoints 2: \t" << kpts2.size() << endl; + cout << "# Matches: \t" << matched1.size() << endl; + cout << "# Inliers: \t" << inliers1.size() << endl; + cout << "# Inliers Ratio: \t" << inlier_ratio << endl; + cout << endl; + + return 0; +}