--- /dev/null
+.. _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<KeyPoint> 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<DMatch> > 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<double>(0) = matched1[i].pt.x;
+ col.at<double>(1) = matched1[i].pt.y;
+
+ col = homography * col;
+ col /= col.at<double>(2);
+ float dist = sqrt( pow(col.at<double>(0) - matched2[i].pt.x, 2) +
+ pow(col.at<double>(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
: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
../feature_flann_matcher/feature_flann_matcher
../feature_homography/feature_homography
../detection_of_planar_objects/detection_of_planar_objects
+ ../akaze_matching/akaze_matching
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.
----
.. 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.
----------
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
-----
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
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
*/
{
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();
CV_PROP bool extended;
CV_PROP bool upright;
+ CV_PROP float threshold;
+ CV_PROP int octaves;
+ CV_PROP int sublevels;
+ CV_PROP int diffusivity;
};
/*!
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();
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 *
*/
#include "precomp.hpp"
-#include "akaze/AKAZEFeatures.h"
+#include "kaze/AKAZEFeatures.h"
+
+#include <iostream>
+using namespace std;
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)
{
}
{
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)
{
{
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:
{
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:
cv::Mat& desc = descriptors.getMatRef();
AKAZEOptions options;
- options.descriptor = static_cast<DESCRIPTOR_TYPE>(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);
img.convertTo(img1_32, CV_32F, 1.0 / 255.0, 0);
AKAZEOptions options;
- options.descriptor = static_cast<DESCRIPTOR_TYPE>(descriptor);
+ options.descriptor = descriptor;
options.descriptor_channels = descriptor_channels;
options.descriptor_size = descriptor_size;
options.img_width = img.cols;
cv::Mat& desc = descriptors.getMatRef();
AKAZEOptions options;
- options.descriptor = static_cast<DESCRIPTOR_TYPE>(descriptor);
+ options.descriptor = descriptor;
options.descriptor_channels = descriptor_channels;
options.descriptor_size = descriptor_size;
options.img_width = img.cols;
CV_Assert((!desc.rows || desc.cols == descriptorSize()));
CV_Assert((!desc.rows || (desc.type() == descriptorType())));
}
-}
\ No newline at end of file
+}
+++ /dev/null
-/**
- * @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<float> 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<cv::KeyPoint>& 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<TEvolution>& ev, const AKAZEOptions& opt)
- : evolution_(&ev)
- , options_(opt)
- {
- }
-
-
- void operator()(const cv::Range& range) const
- {
- std::vector<TEvolution>& 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<TEvolution>* 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<float>(ix)+jx);
- float lxy = *(evolution_[i].Lxy.ptr<float>(ix)+jx);
- float lyy = *(evolution_[i].Lyy.ptr<float>(ix)+jx);
- *(evolution_[i].Ldet.ptr<float>(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<cv::KeyPoint>& 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<cv::KeyPoint> 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<float>(ix)+jx);
-
- // Filter the points with the detector threshold
- if (value > options_.dthreshold && value >= options_.min_dthreshold &&
- value > *(evolution_[i].Ldet.ptr<float>(ix)+jx - 1) &&
- value > *(evolution_[i].Ldet.ptr<float>(ix)+jx + 1) &&
- value > *(evolution_[i].Ldet.ptr<float>(ix - 1) + jx - 1) &&
- value > *(evolution_[i].Ldet.ptr<float>(ix - 1) + jx) &&
- value > *(evolution_[i].Ldet.ptr<float>(ix - 1) + jx + 1) &&
- value > *(evolution_[i].Ldet.ptr<float>(ix + 1) + jx - 1) &&
- value > *(evolution_[i].Ldet.ptr<float>(ix + 1) + jx) &&
- value > *(evolution_[i].Ldet.ptr<float>(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<float>(jx);
- point.pt.y = static_cast<float>(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<cv::KeyPoint>& 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<float>(y)+x + 1)
- - *(evolution_[kpts[i].class_id].Ldet.ptr<float>(y)+x - 1));
- Dy = (0.5f)*(*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y + 1) + x)
- - *(evolution_[kpts[i].class_id].Ldet.ptr<float>(y - 1) + x));
-
- // Compute the Hessian
- Dxx = (*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y)+x + 1)
- + *(evolution_[kpts[i].class_id].Ldet.ptr<float>(y)+x - 1)
- - 2.0f*(*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y)+x)));
-
- Dyy = (*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y + 1) + x)
- + *(evolution_[kpts[i].class_id].Ldet.ptr<float>(y - 1) + x)
- - 2.0f*(*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y)+x)));
-
- Dxy = (0.25f)*(*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y + 1) + x + 1)
- + (*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y - 1) + x - 1)))
- - (0.25f)*(*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y - 1) + x + 1)
- + (*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y + 1) + x - 1)));
-
- // Solve the linear system
- *(A.ptr<float>(0)) = Dxx;
- *(A.ptr<float>(1) + 1) = Dyy;
- *(A.ptr<float>(0) + 1) = *(A.ptr<float>(1)) = Dxy;
- *(b.ptr<float>(0)) = -Dx;
- *(b.ptr<float>(1)) = -Dy;
-
- cv::solve(A, b, dst, DECOMP_LU);
-
- if (fabs(*(dst.ptr<float>(0))) <= 1.0f && fabs(*(dst.ptr<float>(1))) <= 1.0f) {
- kpts[i].pt.x = x + (*(dst.ptr<float>(0)));
- kpts[i].pt.y = y + (*(dst.ptr<float>(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<cv::KeyPoint>& kpts, cv::Mat& desc, std::vector<TEvolution>& 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<float>(i));
- }
- }
-
- void Get_SURF_Descriptor_Upright_64(const cv::KeyPoint& kpt, float* desc) const;
-
-private:
- std::vector<cv::KeyPoint>* keypoints_;
- cv::Mat* descriptors_;
- std::vector<TEvolution>* evolution_;
-};
-
-class SURF_Descriptor_64_Invoker : public cv::ParallelLoopBody
-{
-public:
- SURF_Descriptor_64_Invoker(std::vector<cv::KeyPoint>& kpts, cv::Mat& desc, std::vector<TEvolution>& 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<float>(i));
- }
- }
-
- void Get_SURF_Descriptor_64(const cv::KeyPoint& kpt, float* desc) const;
-
-private:
- std::vector<cv::KeyPoint>* keypoints_;
- cv::Mat* descriptors_;
- std::vector<TEvolution>* evolution_;
-};
-
-class MSURF_Upright_Descriptor_64_Invoker : public cv::ParallelLoopBody
-{
-public:
- MSURF_Upright_Descriptor_64_Invoker(std::vector<cv::KeyPoint>& kpts, cv::Mat& desc, std::vector<TEvolution>& 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<float>(i));
- }
- }
-
- void Get_MSURF_Upright_Descriptor_64(const cv::KeyPoint& kpt, float* desc) const;
-
-private:
- std::vector<cv::KeyPoint>* keypoints_;
- cv::Mat* descriptors_;
- std::vector<TEvolution>* evolution_;
-};
-
-class MSURF_Descriptor_64_Invoker : public cv::ParallelLoopBody
-{
-public:
- MSURF_Descriptor_64_Invoker(std::vector<cv::KeyPoint>& kpts, cv::Mat& desc, std::vector<TEvolution>& 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<float>(i));
- }
- }
-
- void Get_MSURF_Descriptor_64(const cv::KeyPoint& kpt, float* desc) const;
-
-private:
- std::vector<cv::KeyPoint>* keypoints_;
- cv::Mat* descriptors_;
- std::vector<TEvolution>* evolution_;
-};
-
-class Upright_MLDB_Full_Descriptor_Invoker : public cv::ParallelLoopBody
-{
-public:
- Upright_MLDB_Full_Descriptor_Invoker(std::vector<cv::KeyPoint>& kpts, cv::Mat& desc, std::vector<TEvolution>& 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<unsigned char>(i));
- }
- }
-
- void Get_Upright_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned char* desc) const;
-
-private:
- std::vector<cv::KeyPoint>* keypoints_;
- cv::Mat* descriptors_;
- std::vector<TEvolution>* evolution_;
- AKAZEOptions* options_;
-};
-
-class Upright_MLDB_Descriptor_Subset_Invoker : public cv::ParallelLoopBody
-{
-public:
- Upright_MLDB_Descriptor_Subset_Invoker(std::vector<cv::KeyPoint>& kpts,
- cv::Mat& desc,
- std::vector<TEvolution>& 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<unsigned char>(i));
- }
- }
-
- void Get_Upright_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char* desc) const;
-
-private:
- std::vector<cv::KeyPoint>* keypoints_;
- cv::Mat* descriptors_;
- std::vector<TEvolution>* 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<cv::KeyPoint>& kpts, cv::Mat& desc, std::vector<TEvolution>& 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<unsigned char>(i));
- }
- }
-
- void Get_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned char* desc) const;
-
-private:
- std::vector<cv::KeyPoint>* keypoints_;
- cv::Mat* descriptors_;
- std::vector<TEvolution>* evolution_;
- AKAZEOptions* options_;
-};
-
-class MLDB_Descriptor_Subset_Invoker : public cv::ParallelLoopBody
-{
-public:
- MLDB_Descriptor_Subset_Invoker(std::vector<cv::KeyPoint>& kpts,
- cv::Mat& desc,
- std::vector<TEvolution>& 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<unsigned char>(i));
- }
- }
-
- void Get_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char* desc) const;
-
-private:
- std::vector<cv::KeyPoint>* keypoints_;
- cv::Mat* descriptors_;
- std::vector<TEvolution>* 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<cv::KeyPoint>& 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<TEvolution>& 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<float> 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<float>(iy)+ix));
- resY[idx] = gweight*(*(evolution_[level].Ly.ptr<float>(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<TEvolution>& 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<float>(y1)+x1);
- res2 = *(evolution[level].Lx.ptr<float>(y1)+x2);
- res3 = *(evolution[level].Lx.ptr<float>(y2)+x1);
- res4 = *(evolution[level].Lx.ptr<float>(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<float>(y1)+x1);
- res2 = *(evolution[level].Ly.ptr<float>(y1)+x2);
- res3 = *(evolution[level].Ly.ptr<float>(y2)+x1);
- res4 = *(evolution[level].Ly.ptr<float>(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<TEvolution>& 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<float>(y1)+x1);
- res2 = *(evolution[level].Lx.ptr<float>(y1)+x2);
- res3 = *(evolution[level].Lx.ptr<float>(y2)+x1);
- res4 = *(evolution[level].Lx.ptr<float>(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<float>(y1)+x1);
- res2 = *(evolution[level].Ly.ptr<float>(y1)+x2);
- res3 = *(evolution[level].Ly.ptr<float>(y2)+x1);
- res4 = *(evolution[level].Ly.ptr<float>(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<TEvolution>& 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<float>(y1)+x1);
- rx = *(evolution[level].Lx.ptr<float>(y1)+x1);
- ry = *(evolution[level].Ly.ptr<float>(y1)+x1);
-
- di += ri;
- dx += rx;
- dy += ry;
- nsamples++;
- }
- }
-
- di /= nsamples;
- dx /= nsamples;
- dy /= nsamples;
-
- *(values_1.ptr<float>(dcount2)) = di;
- *(values_1.ptr<float>(dcount2)+1) = dx;
- *(values_1.ptr<float>(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<float>(i)) > *(values_1.ptr<float>(j))) {
- desc[dcount1 / 8] |= (1 << (dcount1 % 8));
- }
- dcount1++;
-
- if (*(values_1.ptr<float>(i)+1) > *(values_1.ptr<float>(j)+1)) {
- desc[dcount1 / 8] |= (1 << (dcount1 % 8));
- }
- dcount1++;
-
- if (*(values_1.ptr<float>(i)+2) > *(values_1.ptr<float>(j)+2)) {
- desc[dcount1 / 8] |= (1 << (dcount1 % 8));
- }
- dcount1++;
- }
- }
-
- // Second 3x3 grid
- sample_step = static_cast<int>(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<float>(y1)+x1);
- rx = *(evolution[level].Lx.ptr<float>(y1)+x1);
- ry = *(evolution[level].Ly.ptr<float>(y1)+x1);
-
- di += ri;
- dx += rx;
- dy += ry;
- nsamples++;
- }
- }
-
- di /= nsamples;
- dx /= nsamples;
- dy /= nsamples;
-
- *(values_2.ptr<float>(dcount2)) = di;
- *(values_2.ptr<float>(dcount2)+1) = dx;
- *(values_2.ptr<float>(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<float>(i)) > *(values_2.ptr<float>(j))) {
- desc[dcount1 / 8] |= (1 << (dcount1 % 8));
- }
- dcount1++;
-
- if (*(values_2.ptr<float>(i)+1) > *(values_2.ptr<float>(j)+1)) {
- desc[dcount1 / 8] |= (1 << (dcount1 % 8));
- }
- dcount1++;
-
- if (*(values_2.ptr<float>(i)+2) > *(values_2.ptr<float>(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<float>(y1)+x1);
- rx = *(evolution[level].Lx.ptr<float>(y1)+x1);
- ry = *(evolution[level].Ly.ptr<float>(y1)+x1);
-
- di += ri;
- dx += rx;
- dy += ry;
- nsamples++;
- }
- }
-
- di /= nsamples;
- dx /= nsamples;
- dy /= nsamples;
-
- *(values_3.ptr<float>(dcount2)) = di;
- *(values_3.ptr<float>(dcount2)+1) = dx;
- *(values_3.ptr<float>(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<float>(i)) > *(values_3.ptr<float>(j))) {
- desc[dcount1 / 8] |= (1 << (dcount1 % 8));
- }
- dcount1++;
-
- if (*(values_3.ptr<float>(i)+1) > *(values_3.ptr<float>(j)+1)) {
- desc[dcount1 / 8] |= (1 << (dcount1 % 8));
- }
- dcount1++;
-
- if (*(values_3.ptr<float>(i)+2) > *(values_3.ptr<float>(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<TEvolution>& 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<float>(y1)+x1);
- rx = *(evolution[level].Lx.ptr<float>(y1)+x1);
- ry = *(evolution[level].Ly.ptr<float>(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<float>(dcount2)) = di;
- if (options.descriptor_channels > 1) {
- *(values_1.ptr<float>(dcount2)+1) = dx;
- }
-
- if (options.descriptor_channels > 2) {
- *(values_1.ptr<float>(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<float>(i)) > *(values_1.ptr<float>(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<float>(i)+1) > *(values_1.ptr<float>(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<float>(i)+2) > *(values_1.ptr<float>(j)+2)) {
- desc[dcount1 / 8] |= (1 << (dcount1 % 8));
- }
- dcount1++;
- }
- }
- }
-
- // Second 3x3 grid
- sample_step = static_cast<int>(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<float>(y1)+x1);
- rx = *(evolution[level].Lx.ptr<float>(y1)+x1);
- ry = *(evolution[level].Ly.ptr<float>(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<float>(dcount2)) = di;
- if (options.descriptor_channels > 1) {
- *(values_2.ptr<float>(dcount2)+1) = dx;
- }
-
- if (options.descriptor_channels > 2) {
- *(values_2.ptr<float>(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<float>(i)) > *(values_2.ptr<float>(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<float>(i)+1) > *(values_2.ptr<float>(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<float>(i)+2) > *(values_2.ptr<float>(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<float>(y1)+x1);
- rx = *(evolution[level].Lx.ptr<float>(y1)+x1);
- ry = *(evolution[level].Ly.ptr<float>(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<float>(dcount2)) = di;
- if (options.descriptor_channels > 1)
- *(values_3.ptr<float>(dcount2)+1) = dx;
-
- if (options.descriptor_channels > 2)
- *(values_3.ptr<float>(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<float>(i)) > *(values_3.ptr<float>(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<float>(i)+1) > *(values_3.ptr<float>(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<float>(i)+2) > *(values_3.ptr<float>(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<TEvolution>& 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_<float>::zeros((4 + 9 + 16)*options.descriptor_channels, 1);
-
- // Sample everything, but only do the comparisons
- vector<int> 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<int>(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<float>(y1)+x1);
-
- if (options.descriptor_channels > 1) {
- rx = *(evolution[level].Lx.ptr<float>(y1)+x1);
- ry = *(evolution[level].Ly.ptr<float>(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<float>(options.descriptor_channels*i)) = di;
-
- if (options.descriptor_channels == 2) {
- *(values.ptr<float>(options.descriptor_channels*i + 1)) = dx;
- }
- else if (options.descriptor_channels == 3) {
- *(values.ptr<float>(options.descriptor_channels*i + 1)) = dx;
- *(values.ptr<float>(options.descriptor_channels*i + 2)) = dy;
- }
- }
-
- // Do the comparisons
- const float *vals = values.ptr<float>(0);
- const int *comps = descriptorBits_.ptr<int>(0);
-
- for (int i = 0; i<descriptorBits_.rows; i++) {
- if (vals[comps[2 * 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<TEvolution>& 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_<float>::zeros((4 + 9 + 16)*options.descriptor_channels, 1);
-
- vector<int> steps(3);
- steps.at(0) = options.descriptor_pattern_size;
- steps.at(1) = static_cast<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<int>(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<float>(y1)+x1);
-
- if (options.descriptor_channels > 1) {
- rx = *(evolution[level].Lx.ptr<float>(y1)+x1);
- ry = *(evolution[level].Ly.ptr<float>(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<float>(options.descriptor_channels*i)) = di;
-
- if (options.descriptor_channels == 2) {
- *(values.ptr<float>(options.descriptor_channels*i + 1)) = dx;
- }
- else if (options.descriptor_channels == 3) {
- *(values.ptr<float>(options.descriptor_channels*i + 1)) = dx;
- *(values.ptr<float>(options.descriptor_channels*i + 2)) = dy;
- }
- }
-
- // Do the comparisons
- const float *vals = values.ptr<float>(0);
- const int *comps = descriptorBits_.ptr<int>(0);
-
- for (int i = 0; i<descriptorBits_.rows; i++) {
- if (vals[comps[2 * 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_<int> 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_<int> comps = Mat_<int>(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_<int> samples(29, 3);
- Mat_<int> 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<float>(CV_PI)-atanf(-y / x);
- }
-
- if (x < 0 && y < 0) {
- return static_cast<float>(CV_PI)+atanf(y / x);
- }
-
- if (x >= 0 && y < 0) {
- return static_cast<float>(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
+++ /dev/null
-/**
- * @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<TEvolution> evolution_; ///< Vector of nonlinear diffusion evolution
-
- /// FED parameters
- int ncycles_; ///< Number of cycles
- bool reordering_; ///< Flag for reordering time steps
- std::vector<std::vector<float > > tsteps_; ///< Vector of FED dynamic time steps
- std::vector<int> 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<cv::KeyPoint>& kpts);
- void Compute_Determinant_Hessian_Response(void);
- void Compute_Multiscale_Derivatives(void);
- void Find_Scale_Space_Extrema(std::vector<cv::KeyPoint>& kpts);
- void Do_Subpixel_Refinement(std::vector<cv::KeyPoint>& kpts);
-
- // Feature description methods
- void Compute_Descriptors(std::vector<cv::KeyPoint>& kpts, cv::Mat& desc);
-
- static void Compute_Main_Orientation(cv::KeyPoint& kpt, const std::vector<TEvolution>& 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);
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)
{
}
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);
CV_Assert((!desc.rows || desc.cols == descriptorSize()));
CV_Assert((!desc.rows || (desc.type() == descriptorType())));
}
-}
\ No newline at end of file
+}
* @author Pablo F. Alcantarilla, Jesus Nuevo
*/
-#pragma once
+#ifndef __OPENCV_FEATURES_2D_AKAZE_CONFIG_H__
+#define __OPENCV_FEATURES_2D_AKAZE_CONFIG_H__
/* ************************************************************************* */
// OpenCV
/// 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)
, 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)
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
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
--- /dev/null
+/**
+ * @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 <iostream>
+
+// 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<float> 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<cv::KeyPoint>& 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<TEvolution>& ev, const AKAZEOptions& opt)
+ : evolution_(&ev)
+ , options_(opt)
+ {
+ }
+
+ void operator()(const cv::Range& range) const
+ {
+ std::vector<TEvolution>& 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<TEvolution>* 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<float>(ix)+jx);
+ float lxy = *(evolution_[i].Lxy.ptr<float>(ix)+jx);
+ float lyy = *(evolution_[i].Lyy.ptr<float>(ix)+jx);
+ *(evolution_[i].Ldet.ptr<float>(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<cv::KeyPoint>& 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<cv::KeyPoint> 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<float>(ix)+jx);
+
+ // Filter the points with the detector threshold
+ if (value > options_.dthreshold && value >= options_.min_dthreshold &&
+ value > *(evolution_[i].Ldet.ptr<float>(ix)+jx - 1) &&
+ value > *(evolution_[i].Ldet.ptr<float>(ix)+jx + 1) &&
+ value > *(evolution_[i].Ldet.ptr<float>(ix - 1) + jx - 1) &&
+ value > *(evolution_[i].Ldet.ptr<float>(ix - 1) + jx) &&
+ value > *(evolution_[i].Ldet.ptr<float>(ix - 1) + jx + 1) &&
+ value > *(evolution_[i].Ldet.ptr<float>(ix + 1) + jx - 1) &&
+ value > *(evolution_[i].Ldet.ptr<float>(ix + 1) + jx) &&
+ value > *(evolution_[i].Ldet.ptr<float>(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<float>(jx);
+ point.pt.y = static_cast<float>(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<cv::KeyPoint>& 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<float>(y)+x + 1)
+ - *(evolution_[kpts[i].class_id].Ldet.ptr<float>(y)+x - 1));
+ Dy = (0.5f)*(*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y + 1) + x)
+ - *(evolution_[kpts[i].class_id].Ldet.ptr<float>(y - 1) + x));
+
+ // Compute the Hessian
+ Dxx = (*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y)+x + 1)
+ + *(evolution_[kpts[i].class_id].Ldet.ptr<float>(y)+x - 1)
+ - 2.0f*(*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y)+x)));
+
+ Dyy = (*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y + 1) + x)
+ + *(evolution_[kpts[i].class_id].Ldet.ptr<float>(y - 1) + x)
+ - 2.0f*(*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y)+x)));
+
+ Dxy = (0.25f)*(*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y + 1) + x + 1)
+ + (*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y - 1) + x - 1)))
+ - (0.25f)*(*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y - 1) + x + 1)
+ + (*(evolution_[kpts[i].class_id].Ldet.ptr<float>(y + 1) + x - 1)));
+
+ // Solve the linear system
+ *(A.ptr<float>(0)) = Dxx;
+ *(A.ptr<float>(1) + 1) = Dyy;
+ *(A.ptr<float>(0) + 1) = *(A.ptr<float>(1)) = Dxy;
+ *(b.ptr<float>(0)) = -Dx;
+ *(b.ptr<float>(1)) = -Dy;
+
+ cv::solve(A, b, dst, DECOMP_LU);
+
+ if (fabs(*(dst.ptr<float>(0))) <= 1.0f && fabs(*(dst.ptr<float>(1))) <= 1.0f) {
+ kpts[i].pt.x = x + (*(dst.ptr<float>(0)));
+ kpts[i].pt.y = y + (*(dst.ptr<float>(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<cv::KeyPoint>& kpts, cv::Mat& desc, std::vector<TEvolution>& 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<float>(i));
+ }
+ }
+
+ void Get_SURF_Descriptor_Upright_64(const cv::KeyPoint& kpt, float* desc) const;
+
+private:
+ std::vector<cv::KeyPoint>* keypoints_;
+ cv::Mat* descriptors_;
+ std::vector<TEvolution>* evolution_;
+};
+
+class SURF_Descriptor_64_Invoker : public cv::ParallelLoopBody
+{
+public:
+ SURF_Descriptor_64_Invoker(std::vector<cv::KeyPoint>& kpts, cv::Mat& desc, std::vector<TEvolution>& 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<float>(i));
+ }
+ }
+
+ void Get_SURF_Descriptor_64(const cv::KeyPoint& kpt, float* desc) const;
+
+private:
+ std::vector<cv::KeyPoint>* keypoints_;
+ cv::Mat* descriptors_;
+ std::vector<TEvolution>* evolution_;
+};
+
+class MSURF_Upright_Descriptor_64_Invoker : public cv::ParallelLoopBody
+{
+public:
+ MSURF_Upright_Descriptor_64_Invoker(std::vector<cv::KeyPoint>& kpts, cv::Mat& desc, std::vector<TEvolution>& 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<float>(i));
+ }
+ }
+
+ void Get_MSURF_Upright_Descriptor_64(const cv::KeyPoint& kpt, float* desc) const;
+
+private:
+ std::vector<cv::KeyPoint>* keypoints_;
+ cv::Mat* descriptors_;
+ std::vector<TEvolution>* evolution_;
+};
+
+class MSURF_Descriptor_64_Invoker : public cv::ParallelLoopBody
+{
+public:
+ MSURF_Descriptor_64_Invoker(std::vector<cv::KeyPoint>& kpts, cv::Mat& desc, std::vector<TEvolution>& 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<float>(i));
+ }
+ }
+
+ void Get_MSURF_Descriptor_64(const cv::KeyPoint& kpt, float* desc) const;
+
+private:
+ std::vector<cv::KeyPoint>* keypoints_;
+ cv::Mat* descriptors_;
+ std::vector<TEvolution>* evolution_;
+};
+
+class Upright_MLDB_Full_Descriptor_Invoker : public cv::ParallelLoopBody
+{
+public:
+ Upright_MLDB_Full_Descriptor_Invoker(std::vector<cv::KeyPoint>& kpts, cv::Mat& desc, std::vector<TEvolution>& 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<unsigned char>(i));
+ }
+ }
+
+ void Get_Upright_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned char* desc) const;
+
+private:
+ std::vector<cv::KeyPoint>* keypoints_;
+ cv::Mat* descriptors_;
+ std::vector<TEvolution>* evolution_;
+ AKAZEOptions* options_;
+};
+
+class Upright_MLDB_Descriptor_Subset_Invoker : public cv::ParallelLoopBody
+{
+public:
+ Upright_MLDB_Descriptor_Subset_Invoker(std::vector<cv::KeyPoint>& kpts,
+ cv::Mat& desc,
+ std::vector<TEvolution>& 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<unsigned char>(i));
+ }
+ }
+
+ void Get_Upright_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char* desc) const;
+
+private:
+ std::vector<cv::KeyPoint>* keypoints_;
+ cv::Mat* descriptors_;
+ std::vector<TEvolution>* 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<cv::KeyPoint>& kpts, cv::Mat& desc, std::vector<TEvolution>& 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<unsigned char>(i));
+ }
+ }
+
+ void Get_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned char* desc) const;
+
+private:
+ std::vector<cv::KeyPoint>* keypoints_;
+ cv::Mat* descriptors_;
+ std::vector<TEvolution>* evolution_;
+ AKAZEOptions* options_;
+};
+
+class MLDB_Descriptor_Subset_Invoker : public cv::ParallelLoopBody
+{
+public:
+ MLDB_Descriptor_Subset_Invoker(std::vector<cv::KeyPoint>& kpts,
+ cv::Mat& desc,
+ std::vector<TEvolution>& 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<unsigned char>(i));
+ }
+ }
+
+ void Get_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char* desc) const;
+
+private:
+ std::vector<cv::KeyPoint>* keypoints_;
+ cv::Mat* descriptors_;
+ std::vector<TEvolution>* 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<cv::KeyPoint>& 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<int>(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<TEvolution>& 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<float> 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<float>(iy)+ix));
+ resY[idx] = gweight*(*(evolution_[level].Ly.ptr<float>(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<TEvolution>& 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<float>(y1)+x1);
+ res2 = *(evolution[level].Lx.ptr<float>(y1)+x2);
+ res3 = *(evolution[level].Lx.ptr<float>(y2)+x1);
+ res4 = *(evolution[level].Lx.ptr<float>(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<float>(y1)+x1);
+ res2 = *(evolution[level].Ly.ptr<float>(y1)+x2);
+ res3 = *(evolution[level].Ly.ptr<float>(y2)+x1);
+ res4 = *(evolution[level].Ly.ptr<float>(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<TEvolution>& 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<float>(y1)+x1);
+ res2 = *(evolution[level].Lx.ptr<float>(y1)+x2);
+ res3 = *(evolution[level].Lx.ptr<float>(y2)+x1);
+ res4 = *(evolution[level].Lx.ptr<float>(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<float>(y1)+x1);
+ res2 = *(evolution[level].Ly.ptr<float>(y1)+x2);
+ res3 = *(evolution[level].Ly.ptr<float>(y2)+x1);
+ res4 = *(evolution[level].Ly.ptr<float>(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<TEvolution>& 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<float>(y1)+x1);
+ rx = *(evolution[level].Lx.ptr<float>(y1)+x1);
+ ry = *(evolution[level].Ly.ptr<float>(y1)+x1);
+
+ di += ri;
+ dx += rx;
+ dy += ry;
+ nsamples++;
+ }
+ }
+
+ di /= nsamples;
+ dx /= nsamples;
+ dy /= nsamples;
+
+ *(values_1.ptr<float>(dcount2)) = di;
+ *(values_1.ptr<float>(dcount2)+1) = dx;
+ *(values_1.ptr<float>(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<float>(i)) > *(values_1.ptr<float>(j))) {
+ desc[dcount1 / 8] |= (1 << (dcount1 % 8));
+ }
+ dcount1++;
+
+ if (*(values_1.ptr<float>(i)+1) > *(values_1.ptr<float>(j)+1)) {
+ desc[dcount1 / 8] |= (1 << (dcount1 % 8));
+ }
+ dcount1++;
+
+ if (*(values_1.ptr<float>(i)+2) > *(values_1.ptr<float>(j)+2)) {
+ desc[dcount1 / 8] |= (1 << (dcount1 % 8));
+ }
+ dcount1++;
+ }
+ }
+
+ // Second 3x3 grid
+ sample_step = static_cast<int>(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<float>(y1)+x1);
+ rx = *(evolution[level].Lx.ptr<float>(y1)+x1);
+ ry = *(evolution[level].Ly.ptr<float>(y1)+x1);
+
+ di += ri;
+ dx += rx;
+ dy += ry;
+ nsamples++;
+ }
+ }
+
+ di /= nsamples;
+ dx /= nsamples;
+ dy /= nsamples;
+
+ *(values_2.ptr<float>(dcount2)) = di;
+ *(values_2.ptr<float>(dcount2)+1) = dx;
+ *(values_2.ptr<float>(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<float>(i)) > *(values_2.ptr<float>(j))) {
+ desc[dcount1 / 8] |= (1 << (dcount1 % 8));
+ }
+ dcount1++;
+
+ if (*(values_2.ptr<float>(i)+1) > *(values_2.ptr<float>(j)+1)) {
+ desc[dcount1 / 8] |= (1 << (dcount1 % 8));
+ }
+ dcount1++;
+
+ if (*(values_2.ptr<float>(i)+2) > *(values_2.ptr<float>(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<float>(y1)+x1);
+ rx = *(evolution[level].Lx.ptr<float>(y1)+x1);
+ ry = *(evolution[level].Ly.ptr<float>(y1)+x1);
+
+ di += ri;
+ dx += rx;
+ dy += ry;
+ nsamples++;
+ }
+ }
+
+ di /= nsamples;
+ dx /= nsamples;
+ dy /= nsamples;
+
+ *(values_3.ptr<float>(dcount2)) = di;
+ *(values_3.ptr<float>(dcount2)+1) = dx;
+ *(values_3.ptr<float>(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<float>(i)) > *(values_3.ptr<float>(j))) {
+ desc[dcount1 / 8] |= (1 << (dcount1 % 8));
+ }
+ dcount1++;
+
+ if (*(values_3.ptr<float>(i)+1) > *(values_3.ptr<float>(j)+1)) {
+ desc[dcount1 / 8] |= (1 << (dcount1 % 8));
+ }
+ dcount1++;
+
+ if (*(values_3.ptr<float>(i)+2) > *(values_3.ptr<float>(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<TEvolution>& 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<float>(y1)+x1);
+ rx = *(evolution[level].Lx.ptr<float>(y1)+x1);
+ ry = *(evolution[level].Ly.ptr<float>(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<float>(dcount2)) = di;
+ if (options.descriptor_channels > 1) {
+ *(values_1.ptr<float>(dcount2)+1) = dx;
+ }
+
+ if (options.descriptor_channels > 2) {
+ *(values_1.ptr<float>(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<float>(i)) > *(values_1.ptr<float>(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<float>(i)+1) > *(values_1.ptr<float>(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<float>(i)+2) > *(values_1.ptr<float>(j)+2)) {
+ desc[dcount1 / 8] |= (1 << (dcount1 % 8));
+ }
+ dcount1++;
+ }
+ }
+ }
+
+ // Second 3x3 grid
+ sample_step = static_cast<int>(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<float>(y1)+x1);
+ rx = *(evolution[level].Lx.ptr<float>(y1)+x1);
+ ry = *(evolution[level].Ly.ptr<float>(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<float>(dcount2)) = di;
+ if (options.descriptor_channels > 1) {
+ *(values_2.ptr<float>(dcount2)+1) = dx;
+ }
+
+ if (options.descriptor_channels > 2) {
+ *(values_2.ptr<float>(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<float>(i)) > *(values_2.ptr<float>(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<float>(i)+1) > *(values_2.ptr<float>(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<float>(i)+2) > *(values_2.ptr<float>(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<float>(y1)+x1);
+ rx = *(evolution[level].Lx.ptr<float>(y1)+x1);
+ ry = *(evolution[level].Ly.ptr<float>(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<float>(dcount2)) = di;
+ if (options.descriptor_channels > 1)
+ *(values_3.ptr<float>(dcount2)+1) = dx;
+
+ if (options.descriptor_channels > 2)
+ *(values_3.ptr<float>(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<float>(i)) > *(values_3.ptr<float>(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<float>(i)+1) > *(values_3.ptr<float>(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<float>(i)+2) > *(values_3.ptr<float>(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<TEvolution>& 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_<float>::zeros((4 + 9 + 16)*options.descriptor_channels, 1);
+
+ // Sample everything, but only do the comparisons
+ vector<int> 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<int>(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<float>(y1)+x1);
+
+ if (options.descriptor_channels > 1) {
+ rx = *(evolution[level].Lx.ptr<float>(y1)+x1);
+ ry = *(evolution[level].Ly.ptr<float>(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<float>(options.descriptor_channels*i)) = di;
+
+ if (options.descriptor_channels == 2) {
+ *(values.ptr<float>(options.descriptor_channels*i + 1)) = dx;
+ }
+ else if (options.descriptor_channels == 3) {
+ *(values.ptr<float>(options.descriptor_channels*i + 1)) = dx;
+ *(values.ptr<float>(options.descriptor_channels*i + 2)) = dy;
+ }
+ }
+
+ // Do the comparisons
+ const float *vals = values.ptr<float>(0);
+ const int *comps = descriptorBits_.ptr<int>(0);
+
+ for (int i = 0; i<descriptorBits_.rows; i++) {
+ if (vals[comps[2 * 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<TEvolution>& 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_<float>::zeros((4 + 9 + 16)*options.descriptor_channels, 1);
+
+ vector<int> steps(3);
+ steps.at(0) = options.descriptor_pattern_size;
+ steps.at(1) = static_cast<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<int>(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<float>(y1)+x1);
+
+ if (options.descriptor_channels > 1) {
+ rx = *(evolution[level].Lx.ptr<float>(y1)+x1);
+ ry = *(evolution[level].Ly.ptr<float>(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<float>(options.descriptor_channels*i)) = di;
+
+ if (options.descriptor_channels == 2) {
+ *(values.ptr<float>(options.descriptor_channels*i + 1)) = dx;
+ }
+ else if (options.descriptor_channels == 3) {
+ *(values.ptr<float>(options.descriptor_channels*i + 1)) = dx;
+ *(values.ptr<float>(options.descriptor_channels*i + 2)) = dy;
+ }
+ }
+
+ // Do the comparisons
+ const float *vals = values.ptr<float>(0);
+ const int *comps = descriptorBits_.ptr<int>(0);
+
+ for (int i = 0; i<descriptorBits_.rows; i++) {
+ if (vals[comps[2 * 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_<int> 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_<int> comps = Mat_<int>(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_<int> samples(29, 3);
+ Mat_<int> 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();
+}
--- /dev/null
+/**
+ * @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<TEvolution> evolution_; ///< Vector of nonlinear diffusion evolution
+
+ /// FED parameters
+ int ncycles_; ///< Number of cycles
+ bool reordering_; ///< Flag for reordering time steps
+ std::vector<std::vector<float > > tsteps_; ///< Vector of FED dynamic time steps
+ std::vector<int> 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<cv::KeyPoint>& kpts);
+ void Compute_Determinant_Hessian_Response(void);
+ void Compute_Multiscale_Derivatives(void);
+ void Find_Scale_Space_Extrema(std::vector<cv::KeyPoint>& kpts);
+ void Do_Subpixel_Refinement(std::vector<cv::KeyPoint>& kpts);
+
+ /// Feature description methods
+ void Compute_Descriptors(std::vector<cv::KeyPoint>& kpts, cv::Mat& desc);
+ static void Compute_Main_Orientation(cv::KeyPoint& kpt, const std::vector<TEvolution>& evolution_);
+};
+
+/* ************************************************************************* */
+/// Inline functions
+void generateDescriptorSubsample(cv::Mat& sampleList, cv::Mat& comparisons,
+ int nbits, int pattern_size, int nchannels);
+
+#endif
* @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"
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)
, 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;
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
*/
#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;
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<float> 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<float> 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
// 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
*/
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
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<float>(ix)+jx);
lxy = *(evolution_[i].Lxy.ptr<float>(ix)+jx);
}
}
-//*************************************************************************************
-//*************************************************************************************
-
+/* ************************************************************************* */
/**
* @brief This method selects interesting keypoints through the nonlinear scale space
* @param kpts Vector of keypoints
void KAZEFeatures::Feature_Detection(std::vector<cv::KeyPoint>& 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<TEvolution>& ev) : evolution_(&ev)
+ {
+ }
+
+ void operator()(const cv::Range& range) const
+ {
+ std::vector<TEvolution>& 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<TEvolution>* 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<TEvolution>& ev, std::vector<std::vector<cv::KeyPoint> >& kpts_par,
+ const KAZEOptions& options) : evolution_(&ev), kpts_par_(&kpts_par), options_(options)
+ {
+ }
+
+ void operator()(const cv::Range& range) const
+ {
+ std::vector<TEvolution>& evolution = *evolution_;
+ std::vector<std::vector<cv::KeyPoint> >& 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<float>(ix)+jx);
+
+ // Filter the points with the detector threshold
+ if (value > options_.dthreshold) {
+ if (value >= *(evolution[i].Ldet.ptr<float>(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<float>(evolution[i].sublevel);
+ kpts_par[i - 1].push_back(point);
+ }
+ }
+ }
+ }
+ }
+
+private:
+ std::vector<TEvolution>* evolution_;
+ std::vector<std::vector<cv::KeyPoint> >* 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<cv::KeyPoint>& kpts)
+void KAZEFeatures::Determinant_Hessian(std::vector<cv::KeyPoint>& kpts)
{
int level = 0;
float dist = 0.0, smax = 3.0;
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++) {
}
}
-//*************************************************************************************
-//*************************************************************************************
-
-/**
- * @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<float>(ix)+jx);
-
- // Filter the points with the detector threshold
- if (value > options.dthreshold) {
- if (value >= *(evolution_[level].Ldet.ptr<float>(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
if (fabs(*(dst.ptr<float>(0))) <= 1.0f && fabs(*(dst.ptr<float>(1))) <= 1.0f && fabs(*(dst.ptr<float>(2))) <= 1.0f) {
kpts_[i].pt.x += *(dst.ptr<float>(0));
kpts_[i].pt.y += *(dst.ptr<float>(1));
- dsc = kpts_[i].octave + (kpts_[i].angle + *(dst.ptr<float>(2))) / ((float)(options.nsublevels));
+ dsc = kpts_[i].octave + (kpts_[i].angle + *(dst.ptr<float>(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
}
}
-//*************************************************************************************
-//*************************************************************************************
-
+/* ************************************************************************* */
class KAZE_Descriptor_Invoker : public cv::ParallelLoopBody
{
public:
- KAZE_Descriptor_Invoker(std::vector<cv::KeyPoint> &kpts, cv::Mat &desc, std::vector<TEvolution>& evolution, const KAZEOptions& _options)
- : _kpts(&kpts)
- , _desc(&desc)
- , _evolution(&evolution)
- , options(_options)
+ KAZE_Descriptor_Invoker(std::vector<cv::KeyPoint> &kpts, cv::Mat &desc, std::vector<TEvolution>& evolution, const KAZEOptions& options)
+ : kpts_(&kpts)
+ , desc_(&desc)
+ , evolution_(&evolution)
+ , options_(options)
{
}
void operator() (const cv::Range& range) const
{
- std::vector<cv::KeyPoint> &kpts = *_kpts;
- cv::Mat &desc = *_desc;
- std::vector<TEvolution> &evolution = *_evolution;
+ std::vector<cv::KeyPoint> &kpts = *kpts_;
+ cv::Mat &desc = *desc_;
+ std::vector<TEvolution> &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<float>((int)i));
else
Get_KAZE_Upright_Descriptor_64(kpts[i], desc.ptr<float>((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<float>((int)i));
else
Get_KAZE_Descriptor_64(kpts[i], desc.ptr<float>((int)i));
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<cv::KeyPoint> * _kpts;
- cv::Mat * _desc;
- std::vector<TEvolution> * _evolution;
- KAZEOptions options;
+ std::vector<cv::KeyPoint> * kpts_;
+ cv::Mat * desc_;
+ std::vector<TEvolution> * evolution_;
+ KAZEOptions options_;
};
+/* ************************************************************************* */
/**
* @brief This method computes the set of descriptors through the nonlinear scale space
* @param kpts Vector of keypoints
*/
void KAZEFeatures::Feature_Description(std::vector<cv::KeyPoint> &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<int>(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
}
}
-//*************************************************************************************
-//*************************************************************************************
-
+/* ************************************************************************* */
/**
* @brief This method computes the upright descriptor (not rotation invariant) of
* the provided 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<TEvolution>& evolution_ = *_evolution;
+ std::vector<TEvolution>& evolution = *evolution_;
// Subregion centers for the 4x4 gaussian weighting
float cx = -0.5f, cy = 0.5f;
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<float>(y1)+x1);
- res2 = *(evolution_[level].Lx.ptr<float>(y1)+x2);
- res3 = *(evolution_[level].Lx.ptr<float>(y2)+x1);
- res4 = *(evolution_[level].Lx.ptr<float>(y2)+x2);
+ res1 = *(evolution[level].Lx.ptr<float>(y1)+x1);
+ res2 = *(evolution[level].Lx.ptr<float>(y1)+x2);
+ res3 = *(evolution[level].Lx.ptr<float>(y2)+x1);
+ res4 = *(evolution[level].Lx.ptr<float>(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<float>(y1)+x1);
- res2 = *(evolution_[level].Ly.ptr<float>(y1)+x2);
- res3 = *(evolution_[level].Ly.ptr<float>(y2)+x1);
- res4 = *(evolution_[level].Ly.ptr<float>(y2)+x2);
+ res1 = *(evolution[level].Ly.ptr<float>(y1)+x1);
+ res2 = *(evolution[level].Ly.ptr<float>(y1)+x2);
+ res3 = *(evolution[level].Ly.ptr<float>(y2)+x1);
+ res4 = *(evolution[level].Ly.ptr<float>(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;
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
int kx = 0, ky = 0, i = 0, j = 0, dcount = 0;
int dsize = 0, scale = 0, level = 0;
- std::vector<TEvolution>& evolution_ = *_evolution;
+ std::vector<TEvolution>& evolution = *evolution_;
// Subregion centers for the 4x4 gaussian weighting
float cx = -0.5f, cy = 0.5f;
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<float>(y1)+x1);
- res2 = *(evolution_[level].Lx.ptr<float>(y1)+x2);
- res3 = *(evolution_[level].Lx.ptr<float>(y2)+x1);
- res4 = *(evolution_[level].Lx.ptr<float>(y2)+x2);
+ res1 = *(evolution[level].Lx.ptr<float>(y1)+x1);
+ res2 = *(evolution[level].Lx.ptr<float>(y1)+x2);
+ res3 = *(evolution[level].Lx.ptr<float>(y2)+x1);
+ res4 = *(evolution[level].Lx.ptr<float>(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<float>(y1)+x1);
- res2 = *(evolution_[level].Ly.ptr<float>(y1)+x2);
- res3 = *(evolution_[level].Ly.ptr<float>(y2)+x1);
- res4 = *(evolution_[level].Ly.ptr<float>(y2)+x2);
+ res1 = *(evolution[level].Ly.ptr<float>(y1)+x1);
+ res2 = *(evolution[level].Ly.ptr<float>(y1)+x2);
+ res3 = *(evolution[level].Ly.ptr<float>(y2)+x1);
+ res4 = *(evolution[level].Ly.ptr<float>(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
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
// Subregion centers for the 4x4 gaussian weighting
float cx = -0.5f, cy = 0.5f;
- std::vector<TEvolution>& evolution_ = *_evolution;
+ std::vector<TEvolution>& evolution = *evolution_;
// Set the descriptor size and the sample and pattern sizes
dsize = 128;
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<float>(y1)+x1);
- res2 = *(evolution_[level].Lx.ptr<float>(y1)+x2);
- res3 = *(evolution_[level].Lx.ptr<float>(y2)+x1);
- res4 = *(evolution_[level].Lx.ptr<float>(y2)+x2);
+ res1 = *(evolution[level].Lx.ptr<float>(y1)+x1);
+ res2 = *(evolution[level].Lx.ptr<float>(y1)+x2);
+ res3 = *(evolution[level].Lx.ptr<float>(y2)+x1);
+ res4 = *(evolution[level].Lx.ptr<float>(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<float>(y1)+x1);
- res2 = *(evolution_[level].Ly.ptr<float>(y1)+x2);
- res3 = *(evolution_[level].Ly.ptr<float>(y2)+x1);
- res4 = *(evolution_[level].Ly.ptr<float>(y2)+x2);
+ res1 = *(evolution[level].Ly.ptr<float>(y1)+x1);
+ res2 = *(evolution[level].Ly.ptr<float>(y1)+x2);
+ res3 = *(evolution[level].Ly.ptr<float>(y2)+x1);
+ res4 = *(evolution[level].Ly.ptr<float>(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;
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
int kx = 0, ky = 0, i = 0, j = 0, dcount = 0;
int dsize = 0, scale = 0, level = 0;
- std::vector<TEvolution>& evolution_ = *_evolution;
+ std::vector<TEvolution>& evolution = *evolution_;
// Subregion centers for the 4x4 gaussian weighting
float cx = -0.5f, cy = 0.5f;
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<float>(y1)+x1);
- res2 = *(evolution_[level].Lx.ptr<float>(y1)+x2);
- res3 = *(evolution_[level].Lx.ptr<float>(y2)+x1);
- res4 = *(evolution_[level].Lx.ptr<float>(y2)+x2);
+ res1 = *(evolution[level].Lx.ptr<float>(y1)+x1);
+ res2 = *(evolution[level].Lx.ptr<float>(y1)+x2);
+ res3 = *(evolution[level].Lx.ptr<float>(y2)+x1);
+ res4 = *(evolution[level].Lx.ptr<float>(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<float>(y1)+x1);
- res2 = *(evolution_[level].Ly.ptr<float>(y1)+x2);
- res3 = *(evolution_[level].Ly.ptr<float>(y2)+x1);
- res4 = *(evolution_[level].Ly.ptr<float>(y2)+x2);
+ res1 = *(evolution[level].Ly.ptr<float>(y1)+x1);
+ res2 = *(evolution[level].Ly.ptr<float>(y1)+x2);
+ res3 = *(evolution[level].Ly.ptr<float>(y2)+x1);
+ res4 = *(evolution[level].Ly.ptr<float>(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
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<float>(i)+j) = *(c.ptr<float>(i)+j) + *(c.ptr<float>(i + 1) + j);
- }
- }
-
- for (int j = 0; j < py_.cols; j++) {
- *(py_.ptr<float>(0) + j) = *(qr_.ptr<float>(0) + j);
- }
-
- for (int j = 0; j < py_.cols; j++) {
- *(py_.ptr<float>(py_.rows - 1) + j) = *(qr_.ptr<float>(qr_.rows - 1) + j);
- }
-
- for (int i = 1; i < py_.rows - 1; i++) {
- for (int j = 0; j < py_.cols; j++) {
- *(py_.ptr<float>(i)+j) = *(qr_.ptr<float>(i - 1) + j) + *(qr_.ptr<float>(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<float>(i)+j) = *(c.ptr<float>(i)+j) + *(c.ptr<float>(i)+j + 1);
- }
- }
-
- for (int i = 0; i < px_.rows; i++) {
- *(px_.ptr<float>(i)) = *(qc_.ptr<float>(i));
- }
-
- for (int i = 0; i < px_.rows; i++) {
- *(px_.ptr<float>(i)+px_.cols - 1) = *(qc_.ptr<float>(i)+qc_.cols - 1);
- }
-
- for (int j = 1; j < px_.cols - 1; j++) {
- for (int i = 0; i < px_.rows; i++) {
- *(px_.ptr<float>(i)+j) = *(qc_.ptr<float>(i)+j - 1) + *(qc_.ptr<float>(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<float>(0) + j) = *(a.ptr<float>(0) + j);
- }
-
- for (int j = 0; j < y.cols; j++) {
- *(y.ptr<float>(0) + j) = *(Ld.ptr<float>(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<float>(k - 1) + j) = *(b.ptr<float>(k - 1) + j) / *(m.ptr<float>(k - 1) + j);
- }
-
- for (int j = 0; j < m.cols; j++) {
- *(m.ptr<float>(k)+j) = *(a.ptr<float>(k)+j) - *(l.ptr<float>(k - 1) + j)*(*(b.ptr<float>(k - 1) + j));
- }
-
- for (int j = 0; j < y.cols; j++) {
- *(y.ptr<float>(k)+j) = *(Ld.ptr<float>(k)+j) - *(l.ptr<float>(k - 1) + j)*(*(y.ptr<float>(k - 1) + j));
- }
- }
-
- // 2. Backward substitution U*x = y
- for (int j = 0; j < y.cols; j++) {
- *(x.ptr<float>(n - 1) + j) = (*(y.ptr<float>(n - 1) + j)) / (*(m.ptr<float>(n - 1) + j));
- }
-
- for (int i = n - 2; i >= 0; i--) {
- for (int j = 0; j < x.cols; j++) {
- *(x.ptr<float>(i)+j) = (*(y.ptr<float>(i)+j) - (*(b.ptr<float>(i)+j))*(*(x.ptr<float>(i + 1) + j))) / (*(m.ptr<float>(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<float>(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);
}
* @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<TEvolution> evolution_; // Vector of nonlinear diffusion evolution
+ /// Parameters of the Nonlinear diffusion class
+ KAZEOptions options_; ///< Configuration options for KAZE
+ std::vector<TEvolution> 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<std::vector<cv::KeyPoint> > kpts_par_;
- // FED parameters
- int ncycles_; // Number of cycles
- bool reordering_; // Flag for reordering time steps
- std::vector<std::vector<float > > tsteps_; // Vector of FED dynamic time steps
- std::vector<int> 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<std::vector<float > > tsteps_; ///< Vector of FED dynamic time steps
+ std::vector<int> 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<cv::KeyPoint>& kpts);
void Feature_Description(std::vector<cv::KeyPoint>& kpts, cv::Mat& desc);
-
static void Compute_Main_Orientation(cv::KeyPoint& kpt, const std::vector<TEvolution>& 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<cv::KeyPoint>& kpts);
- void Find_Extremum_Threading(const int& level);
+ void Determinant_Hessian(std::vector<cv::KeyPoint>& kpts);
void Do_Subpixel_Refinement(std::vector<cv::KeyPoint>& 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
--- /dev/null
+/**
+ * @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
-#ifndef FED_H
-#define FED_H
+#ifndef __OPENCV_FEATURES_2D_FED_H__
+#define __OPENCV_FEATURES_2D_FED_H__
//******************************************************************************
//******************************************************************************
//*************************************************************************************
//*************************************************************************************
-#endif // FED_H
+#endif // __OPENCV_FEATURES_2D_FED_H__
* @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
--- /dev/null
+#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<float>(CV_PI)-atanf(-y / x);
+ }
+
+ if (x < 0 && y < 0) {
+ return static_cast<float>(CV_PI)+atanf(y / x);
+ }
+
+ if (x >= 0 && y < 0) {
+ return static_cast<float>(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
typedef typename Distance::ResultType DistanceType;
CV_DescriptorExtractorTest( const string _name, DistanceType _maxDist, const Ptr<DescriptorExtractor>& _dextractor,
- Distance d = Distance() ):
- name(_name), maxDist(_maxDist), dextractor(_dextractor), distance(d) {}
+ Distance d = Distance(), Ptr<FeatureDetector> _detector = Ptr<FeatureDetector>()):
+ name(_name), maxDist(_maxDist), dextractor(_dextractor), distance(d) , detector(_detector) {}
+
+ ~CV_DescriptorExtractorTest()
+ {
+ if(!detector.empty())
+ detector.release();
+ }
protected:
virtual void createDescriptorExtractor() {}
// Read the test image.
string imgFilename = string(ts->get_data_path()) + FEATURES2D_DIR + "/" + IMAGE_FILENAME;
-
Mat img = imread( imgFilename );
if( img.empty() )
{
ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_TEST_DATA );
return;
}
-
vector<KeyPoint> 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 );
}
}
}
- 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 );
const DistanceType maxDist;
Ptr<DescriptorExtractor> dextractor;
Distance distance;
+ Ptr<FeatureDetector> detector;
private:
CV_DescriptorExtractorTest& operator=(const CV_DescriptorExtractorTest&) { return *this; }
DescriptorExtractor::create("OpponentBRIEF") );
test.safe_run();
}
+
+TEST( Features2d_DescriptorExtractor_KAZE, regression )
+{
+ CV_DescriptorExtractorTest< L2<float> > test( "descriptor-kaze", 0.03f,
+ DescriptorExtractor::create("KAZE"),
+ L2<float>(), FeatureDetector::create("KAZE"));
+ test.safe_run();
+}
+
+TEST( Features2d_DescriptorExtractor_AKAZE, regression )
+{
+ CV_DescriptorExtractorTest<Hamming> test( "descriptor-akaze", (CV_DescriptorExtractorTest<Hamming>::DistanceType)12.f,
+ DescriptorExtractor::create("AKAZE"),
+ Hamming(), FeatureDetector::create("AKAZE"));
+ test.safe_run();
+}
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") );
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<FeatureDetector>("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<FeatureDetector>(new cv::AKAZE(cv::AKAZE::DESCRIPTOR_KAZE)));
+ CV_FeatureDetectorKeypointsTest test_kaze(cv::Ptr<FeatureDetector>(new cv::AKAZE(cv::DESCRIPTOR_KAZE)));
test_kaze.safe_run();
- CV_FeatureDetectorKeypointsTest test_mldb(cv::Ptr<FeatureDetector>(new cv::AKAZE(cv::AKAZE::DESCRIPTOR_MLDB)));
+ CV_FeatureDetectorKeypointsTest test_mldb(cv::Ptr<FeatureDetector>(new cv::AKAZE(cv::DESCRIPTOR_MLDB)));
test_mldb.safe_run();
}
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<FeatureDetector>("Feature2D.KAZE"),
0.08f,
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<FeatureDetector>("Feature2D.AKAZE"),
0.08f,
--- /dev/null
+<?xml version="1.0"?>
+<opencv_storage>
+<H13 type_id="opencv-matrix">
+ <rows>3</rows>
+ <cols>3</cols>
+ <dt>d</dt>
+ <data>
+ 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 </data></H13>
+</opencv_storage>
--- /dev/null
+#include <opencv2/features2d.hpp>
+#include <opencv2/imgcodecs.hpp>
+#include <opencv2/opencv.hpp>
+#include <vector>
+#include <iostream>
+
+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<KeyPoint> kpts1, kpts2;
+ Mat desc1, desc2;
+
+ AKAZE akaze;
+ akaze(img1, noArray(), kpts1, desc1);
+ akaze(img2, noArray(), kpts2, desc2);
+
+ BFMatcher matcher(NORM_HAMMING);
+ vector< vector<DMatch> > nn_matches;
+ matcher.knnMatch(desc1, desc2, nn_matches, 2);
+
+ vector<KeyPoint> matched1, matched2, inliers1, inliers2;
+ vector<DMatch> 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<double>(0) = matched1[i].pt.x;
+ col.at<double>(1) = matched1[i].pt.y;
+
+ col = homography * col;
+ col /= col.at<double>(2);
+ double dist = sqrt( pow(col.at<double>(0) - matched2[i].pt.x, 2) +
+ pow(col.at<double>(1) - matched2[i].pt.y, 2));
+
+ if(dist < inlier_threshold) {
+ int new_i = static_cast<int>(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;
+}
projects / profile / ivi / opencv.git / commitdiff
