Skip to content

O
opencv
Commit 7f829608 authored by f-morozov's avatar f-morozov
Browse files
Options

AKAZE fixes, tests and tutorial

parent 27780248
Hide whitespace changes
Inline Side-by-side
Showing with 2755 additions and 2771 deletions
doc/tutorials/features2d/akaze_matching/akaze_matching.rst 0 → 100644
View file @ 7f829608
.. _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
doc/tutorials/features2d/table_of_content_features2d/table_of_content_features2d.rst
View file @ 7f829608
......@@ -183,6 +183,25 @@ Learn about how to use the feature points detectors, descriptors and matching f
:height: 90pt
:width: 90pt
+
.. tabularcolumns:: m{100pt} m{300pt}
.. cssclass:: toctableopencv
===================== ==============================================
|AkazeMatch| **Title:** :ref:`akazeMatching`
*Compatibility:* > OpenCV 3.0
*Author:* Fedor Morozov
Use *AKAZE* local features to find correspondence between two images.
===================== ==============================================
.. |AkazeMatch| image:: images/AKAZE_Match_Tutorial_Cover.png
:height: 90pt
:width: 90pt
.. raw:: latex
\pagebreak
......@@ -201,3 +220,4 @@ Learn about how to use the feature points detectors, descriptors and matching f
../feature_flann_matcher/feature_flann_matcher
../feature_homography/feature_homography
../detection_of_planar_objects/detection_of_planar_objects
../akaze_matching/akaze_matching
modules/features2d/doc/feature_detection_and_description.rst
View file @ 7f829608
......@@ -31,7 +31,7 @@ Detects corners using the FAST algorithm
Detects corners using the FAST algorithm by [Rosten06]_.
..note:: In Python API, types are given as ``cv2.FAST_FEATURE_DETECTOR_TYPE_5_8``, ``cv2.FAST_FEATURE_DETECTOR_TYPE_7_12`` and ``cv2.FAST_FEATURE_DETECTOR_TYPE_9_16``. For corner detection, use ``cv2.FAST.detect()`` method.
.. note:: In Python API, types are given as ``cv2.FAST_FEATURE_DETECTOR_TYPE_5_8``, ``cv2.FAST_FEATURE_DETECTOR_TYPE_7_12`` and ``cv2.FAST_FEATURE_DETECTOR_TYPE_9_16``. For corner detection, use ``cv2.FAST.detect()`` method.
.. [Rosten06] E. Rosten. Machine Learning for High-speed Corner Detection, 2006.
......@@ -254,7 +254,17 @@ KAZE
----
.. ocv:class:: KAZE : public Feature2D
Class implementing the KAZE keypoint detector and descriptor extractor, described in [ABD12]_.
Class implementing the KAZE keypoint detector and descriptor extractor, described in [ABD12]_. ::
class CV_EXPORTS_W KAZE : public Feature2D
{
public:
CV_WRAP KAZE();
CV_WRAP explicit KAZE(bool extended, bool upright, float threshold = 0.001f,
int octaves = 4, int sublevels = 4, int diffusivity = DIFF_PM_G2);
};
.. note:: AKAZE descriptor can only be used with KAZE or AKAZE keypoints
.. [ABD12] KAZE Features. Pablo F. Alcantarilla, Adrien Bartoli and Andrew J. Davison. In European Conference on Computer Vision (ECCV), Fiorenze, Italy, October 2012.
......@@ -262,12 +272,14 @@ KAZE::KAZE
----------
The KAZE constructor
.. ocv:function:: KAZE::KAZE(bool extended, bool upright)
.. ocv:function:: KAZE::KAZE(bool extended, bool upright, float threshold, int octaves, int sublevels, int diffusivity)
:param extended: Set to enable extraction of extended (128-byte) descriptor.
:param upright: Set to enable use of upright descriptors (non rotation-invariant).
:param threshold: Detector response threshold to accept point
:param octaves: Maximum octave evolution of the image
:param sublevels: Default number of sublevels per scale level
:param diffusivity: Diffusivity type. DIFF_PM_G1, DIFF_PM_G2, DIFF_WEICKERT or DIFF_CHARBONNIER
AKAZE
-----
......@@ -278,25 +290,25 @@ Class implementing the AKAZE keypoint detector and descriptor extractor, describ
class CV_EXPORTS_W AKAZE : public Feature2D
{
public:
/// AKAZE Descriptor Type
enum DESCRIPTOR_TYPE {
DESCRIPTOR_KAZE_UPRIGHT = 2, ///< Upright descriptors, not invariant to rotation
DESCRIPTOR_KAZE = 3,
DESCRIPTOR_MLDB_UPRIGHT = 4, ///< Upright descriptors, not invariant to rotation
DESCRIPTOR_MLDB = 5
};
CV_WRAP AKAZE();
explicit AKAZE(DESCRIPTOR_TYPE descriptor_type, int descriptor_size = 0, int descriptor_channels = 3);
CV_WRAP explicit AKAZE(int descriptor_type, int descriptor_size = 0, int descriptor_channels = 3,
float threshold = 0.001f, int octaves = 4, int sublevels = 4, int diffusivity = DIFF_PM_G2);
};
.. note:: AKAZE descriptor can only be used with KAZE or AKAZE keypoints
.. [ANB13] Fast Explicit Diffusion for Accelerated Features in Nonlinear Scale Spaces. Pablo F. Alcantarilla, Jesús Nuevo and Adrien Bartoli. In British Machine Vision Conference (BMVC), Bristol, UK, September 2013.
AKAZE::AKAZE
------------
The AKAZE constructor
.. ocv:function:: AKAZE::AKAZE(DESCRIPTOR_TYPE descriptor_type, int descriptor_size = 0, int descriptor_channels = 3)
.. ocv:function:: AKAZE::AKAZE(int descriptor_type, int descriptor_size, int descriptor_channels, float threshold, int octaves, int sublevels, int diffusivity)
:param descriptor_type: Type of the extracted descriptor.
:param descriptor_type: Type of the extracted descriptor: DESCRIPTOR_KAZE, DESCRIPTOR_KAZE_UPRIGHT, DESCRIPTOR_MLDB or DESCRIPTOR_MLDB_UPRIGHT.
:param descriptor_size: Size of the descriptor in bits. 0 -> Full size
:param descriptor_channels: Number of channels in the descriptor (1, 2, 3).
:param descriptor_channels: Number of channels in the descriptor (1, 2, 3)
:param threshold: Detector response threshold to accept point
:param octaves: Maximum octave evolution of the image
:param sublevels: Default number of sublevels per scale level
:param diffusivity: Diffusivity type. DIFF_PM_G1, DIFF_PM_G2, DIFF_WEICKERT or DIFF_CHARBONNIER
modules/features2d/include/opencv2/features2d.hpp
View file @ 7f829608
......@@ -895,6 +895,22 @@ protected:
PixelTestFn test_fn_;
};
// KAZE/AKAZE diffusivity
enum {
DIFF_PM_G1 = 0,
DIFF_PM_G2 = 1,
DIFF_WEICKERT = 2,
DIFF_CHARBONNIER = 3
};
// AKAZE descriptor type
enum {
DESCRIPTOR_KAZE_UPRIGHT = 2, ///< Upright descriptors, not invariant to rotation
DESCRIPTOR_KAZE = 3,
DESCRIPTOR_MLDB_UPRIGHT = 4, ///< Upright descriptors, not invariant to rotation
DESCRIPTOR_MLDB = 5
};
/*!
KAZE implementation
*/
......@@ -902,7 +918,8 @@ class CV_EXPORTS_W KAZE : public Feature2D
{
public:
CV_WRAP KAZE();
CV_WRAP explicit KAZE(bool extended, bool upright);
CV_WRAP explicit KAZE(bool extended, bool upright, float threshold = 0.001f,
int octaves = 4, int sublevels = 4, int diffusivity = DIFF_PM_G2);
virtual ~KAZE();
......@@ -928,6 +945,10 @@ protected:
CV_PROP bool extended;
CV_PROP bool upright;
CV_PROP float threshold;
CV_PROP int octaves;
CV_PROP int sublevels;
CV_PROP int diffusivity;
};
/*!
......@@ -936,16 +957,9 @@ AKAZE implementation
class CV_EXPORTS_W AKAZE : public Feature2D
{
public:
/// AKAZE Descriptor Type
enum DESCRIPTOR_TYPE {
DESCRIPTOR_KAZE_UPRIGHT = 2, ///< Upright descriptors, not invariant to rotation
DESCRIPTOR_KAZE = 3,
DESCRIPTOR_MLDB_UPRIGHT = 4, ///< Upright descriptors, not invariant to rotation
DESCRIPTOR_MLDB = 5
};
CV_WRAP AKAZE();
explicit AKAZE(DESCRIPTOR_TYPE descriptor_type, int descriptor_size = 0, int descriptor_channels = 3);
CV_WRAP explicit AKAZE(int descriptor_type, int descriptor_size = 0, int descriptor_channels = 3,
float threshold = 0.001f, int octaves = 4, int sublevels = 4, int diffusivity = DIFF_PM_G2);
virtual ~AKAZE();
......@@ -973,7 +987,10 @@ protected:
CV_PROP int descriptor;
CV_PROP int descriptor_channels;
CV_PROP int descriptor_size;
CV_PROP float threshold;
CV_PROP int octaves;
CV_PROP int sublevels;
CV_PROP int diffusivity;
};
/****************************************************************************************\
* Distance *
......
modules/features2d/src/akaze.cpp
View file @ 7f829608
......@@ -49,7 +49,10 @@ http://www.robesafe.com/personal/pablo.alcantarilla/papers/Alcantarilla13bmvc.pd
*/
#include "precomp.hpp"
#include "akaze/AKAZEFeatures.h"
#include "kaze/AKAZEFeatures.h"
#include <iostream>
using namespace std;
namespace cv
{
......@@ -57,13 +60,22 @@ namespace cv
: descriptor(DESCRIPTOR_MLDB)
, descriptor_channels(3)
, descriptor_size(0)
, threshold(0.001f)
, octaves(4)
, sublevels(4)
, diffusivity(DIFF_PM_G2)
{
}
AKAZE::AKAZE(DESCRIPTOR_TYPE _descriptor_type, int _descriptor_size, int _descriptor_channels)
AKAZE::AKAZE(int _descriptor_type, int _descriptor_size, int _descriptor_channels,
float _threshold, int _octaves, int _sublevels, int _diffusivity)
: descriptor(_descriptor_type)
, descriptor_channels(_descriptor_channels)
, descriptor_size(_descriptor_size)
, threshold(_threshold)
, octaves(_octaves)
, sublevels(_sublevels)
, diffusivity(_diffusivity)
{
}
......@@ -78,12 +90,12 @@ namespace cv
{
switch (descriptor)
{
case cv::AKAZE::DESCRIPTOR_KAZE:
case cv::AKAZE::DESCRIPTOR_KAZE_UPRIGHT:
case cv::DESCRIPTOR_KAZE:
case cv::DESCRIPTOR_KAZE_UPRIGHT:
return 64;
case cv::AKAZE::DESCRIPTOR_MLDB:
case cv::AKAZE::DESCRIPTOR_MLDB_UPRIGHT:
case cv::DESCRIPTOR_MLDB:
case cv::DESCRIPTOR_MLDB_UPRIGHT:
// We use the full length binary descriptor -> 486 bits
if (descriptor_size == 0)
{
......@@ -106,12 +118,12 @@ namespace cv
{
switch (descriptor)
{
case cv::AKAZE::DESCRIPTOR_KAZE:
case cv::AKAZE::DESCRIPTOR_KAZE_UPRIGHT:
case cv::DESCRIPTOR_KAZE:
case cv::DESCRIPTOR_KAZE_UPRIGHT:
return CV_32F;
case cv::AKAZE::DESCRIPTOR_MLDB:
case cv::AKAZE::DESCRIPTOR_MLDB_UPRIGHT:
case cv::DESCRIPTOR_MLDB:
case cv::DESCRIPTOR_MLDB_UPRIGHT:
return CV_8U;
default:
......@@ -124,12 +136,12 @@ namespace cv
{
switch (descriptor)
{
case cv::AKAZE::DESCRIPTOR_KAZE:
case cv::AKAZE::DESCRIPTOR_KAZE_UPRIGHT:
case cv::DESCRIPTOR_KAZE:
case cv::DESCRIPTOR_KAZE_UPRIGHT:
return cv::NORM_L2;
case cv::AKAZE::DESCRIPTOR_MLDB:
case cv::AKAZE::DESCRIPTOR_MLDB_UPRIGHT:
case cv::DESCRIPTOR_MLDB:
case cv::DESCRIPTOR_MLDB_UPRIGHT:
return cv::NORM_HAMMING;
default:
......@@ -153,11 +165,15 @@ namespace cv
cv::Mat& desc = descriptors.getMatRef();
AKAZEOptions options;
options.descriptor = static_cast<DESCRIPTOR_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);
......@@ -188,7 +204,7 @@ namespace cv
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;
......@@ -216,7 +232,7 @@ namespace cv
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;
......@@ -229,4 +245,4 @@ namespace cv
CV_Assert((!desc.rows || desc.cols == descriptorSize()));
CV_Assert((!desc.rows || (desc.type() == descriptorType())));
}
}
\ No newline at end of file
}
modules/features2d/src/akaze/AKAZEFeatures.cpp deleted 100644 → 0
View file @ 27780248
/**
* @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
modules/features2d/src/akaze/AKAZEFeatures.h deleted 100644 → 0
View file @ 27780248
/**
* @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);
modules/features2d/src/kaze.cpp
View file @ 7f829608
......@@ -55,12 +55,21 @@ namespace cv
KAZE::KAZE()
: extended(false)
, upright(false)
, threshold(0.001f)
, octaves(4)
, sublevels(4)
, diffusivity(DIFF_PM_G2)
{
}
KAZE::KAZE(bool _extended, bool _upright)
KAZE::KAZE(bool _extended, bool _upright, float _threshold, int _octaves,
int _sublevels, int _diffusivity)
: extended(_extended)
, upright(_upright)
, threshold(_threshold)
, octaves(_octaves)
, sublevels(_sublevels)
, diffusivity(_diffusivity)
{
}
......@@ -111,6 +120,10 @@ namespace cv
options.img_height = img.rows;
options.extended = extended;
options.upright = upright;
options.dthreshold = threshold;
options.omax = octaves;
options.nsublevels = sublevels;
options.diffusivity = diffusivity;
KAZEFeatures impl(options);
impl.Create_Nonlinear_Scale_Space(img1_32);
......@@ -180,4 +193,4 @@ namespace cv
CV_Assert((!desc.rows || desc.cols == descriptorSize()));
CV_Assert((!desc.rows || (desc.type() == descriptorType())));
}
}
\ No newline at end of file
}
modules/features2d/src/akaze/AKAZEConfig.h → modules/features2d/src/kaze/AKAZEConfig.h
View file @ 7f829608
......@@ -5,7 +5,8 @@
* @author Pablo F. Alcantarilla, Jesus Nuevo
*/
#pragma once
#ifndef __OPENCV_FEATURES_2D_AKAZE_CONFIG_H__
#define __OPENCV_FEATURES_2D_AKAZE_CONFIG_H__
/* ************************************************************************* */
// OpenCV
......@@ -28,14 +29,6 @@ const float gauss25[7][7] = {
/// AKAZE configuration options structure
struct AKAZEOptions {
/// AKAZE Diffusivities
enum DIFFUSIVITY_TYPE {
PM_G1 = 0,
PM_G2 = 1,
WEICKERT = 2,
CHARBONNIER = 3
};
AKAZEOptions()
: omax(4)
, nsublevels(4)
......@@ -44,12 +37,12 @@ struct AKAZEOptions {
, soffset(1.6f)
, derivative_factor(1.5f)
, sderivatives(1.0)
, diffusivity(PM_G2)
, diffusivity(cv::DIFF_PM_G2)
, dthreshold(0.001f)
, min_dthreshold(0.00001f)
, descriptor(cv::AKAZE::DESCRIPTOR_MLDB)
, descriptor(cv::DESCRIPTOR_MLDB)
, descriptor_size(0)
, descriptor_channels(3)
, descriptor_pattern_size(10)
......@@ -67,12 +60,12 @@ struct AKAZEOptions {
float soffset; ///< Base scale offset (sigma units)
float derivative_factor; ///< Factor for the multiscale derivatives
float sderivatives; ///< Smoothing factor for the derivatives
DIFFUSIVITY_TYPE diffusivity; ///< Diffusivity type
int diffusivity; ///< Diffusivity type
float dthreshold; ///< Detector response threshold to accept point
float min_dthreshold; ///< Minimum detector threshold to accept a point
cv::AKAZE::DESCRIPTOR_TYPE descriptor; ///< Type of descriptor
int descriptor; ///< Type of descriptor
int descriptor_size; ///< Size of the descriptor in bits. 0->Full size
int descriptor_channels; ///< Number of channels in the descriptor (1, 2, 3)
int descriptor_pattern_size; ///< Actual patch size is 2*pattern_size*point.scale
......@@ -82,28 +75,4 @@ struct AKAZEOptions {
int kcontrast_nbins; ///< Number of bins for the contrast factor histogram
};
/* ************************************************************************* */
/// AKAZE nonlinear diffusion filtering evolution
struct TEvolution {
TEvolution() {
etime = 0.0f;
esigma = 0.0f;
octave = 0;
sublevel = 0;
sigma_size = 0;
}
cv::Mat Lx, Ly; // First order spatial derivatives
cv::Mat Lxx, Lxy, Lyy; // Second order spatial derivatives
cv::Mat Lflow; // Diffusivity image
cv::Mat Lt; // Evolution image
cv::Mat Lsmooth; // Smoothed image
cv::Mat Lstep; // Evolution step update
cv::Mat Ldet; // Detector response
float etime; // Evolution time
float esigma; // Evolution sigma. For linear diffusion t = sigma^2 / 2
size_t octave; // Image octave
size_t sublevel; // Image sublevel in each octave
size_t sigma_size; // Integer sigma. For computing the feature detector responses
};
\ No newline at end of file
#endif
modules/features2d/src/kaze/AKAZEFeatures.cpp 0 → 100644
View file @ 7f829608
/**
* @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();
}
modules/features2d/src/kaze/AKAZEFeatures.h 0 → 100644
View file @ 7f829608
/**
* @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
modules/features2d/src/kaze/KAZEConfig.h
View file @ 7f829608
......@@ -5,7 +5,8 @@
* @author Pablo F. Alcantarilla
*/
#pragma once
#ifndef __OPENCV_FEATURES_2D_AKAZE_CONFIG_H__
#define __OPENCV_FEATURES_2D_AKAZE_CONFIG_H__
// OpenCV Includes
#include "precomp.hpp"
......@@ -15,14 +16,8 @@
struct KAZEOptions {
enum DIFFUSIVITY_TYPE {
PM_G1 = 0,
PM_G2 = 1,
WEICKERT = 2
};
KAZEOptions()
: diffusivity(PM_G2)
: diffusivity(cv::DIFF_PM_G2)
, soffset(1.60f)
, omax(4)
......@@ -33,20 +28,13 @@ struct KAZEOptions {
, dthreshold(0.001f)
, kcontrast(0.01f)
, kcontrast_percentille(0.7f)
, kcontrast_bins(300)
, use_fed(true)
, kcontrast_bins(300)
, upright(false)
, extended(false)
, use_clipping_normalilzation(false)
, clipping_normalization_ratio(1.6f)
, clipping_normalization_niter(5)
{
}
DIFFUSIVITY_TYPE diffusivity;
int diffusivity;
float soffset;
int omax;
int nsublevels;
......@@ -57,27 +45,8 @@ struct KAZEOptions {
float kcontrast;
float kcontrast_percentille;
int kcontrast_bins;
bool use_fed;
bool upright;
bool extended;
bool use_clipping_normalilzation;
float clipping_normalization_ratio;
int clipping_normalization_niter;
};
struct TEvolution {
cv::Mat Lx, Ly; // First order spatial derivatives
cv::Mat Lxx, Lxy, Lyy; // Second order spatial derivatives
cv::Mat Lflow; // Diffusivity image
cv::Mat Lt; // Evolution image
cv::Mat Lsmooth; // Smoothed image
cv::Mat Lstep; // Evolution step update
cv::Mat Ldet; // Detector response
float etime; // Evolution time
float esigma; // Evolution sigma. For linear diffusion t = sigma^2 / 2
float octave; // Image octave
float sublevel; // Image sublevel in each octave
int sigma_size; // Integer esigma. For computing the feature detector responses
};
#endif
modules/features2d/src/kaze/KAZEFeatures.cpp
View file @ 7f829608
......@@ -22,22 +22,21 @@
*/
#include "KAZEFeatures.h"
#include "utils.h"
// Namespaces
using namespace std;
using namespace cv;
using namespace cv::details::kaze;
//*******************************************************************************
//*******************************************************************************
/* ************************************************************************* */
/**
* @brief KAZE constructor with input options
* @param options KAZE configuration options
* @note The constructor allocates memory for the nonlinear scale space
*/
KAZEFeatures::KAZEFeatures(KAZEOptions& _options)
: options(_options)
KAZEFeatures::KAZEFeatures(KAZEOptions& options)
: options_(options)
{
ncycles_ = 0;
reordering_ = true;
......@@ -46,70 +45,48 @@ KAZEFeatures::KAZEFeatures(KAZEOptions& _options)
Allocate_Memory_Evolution();
}
//*******************************************************************************
//*******************************************************************************
/* ************************************************************************* */
/**
* @brief This method allocates the memory for the nonlinear diffusion evolution
*/
void KAZEFeatures::Allocate_Memory_Evolution(void) {
// Allocate the dimension of the matrices for the evolution
for (int i = 0; i <= options.omax - 1; i++) {
for (int j = 0; j <= options.nsublevels - 1; j++) {
for (int i = 0; i <= options_.omax - 1; i++) {
for (int j = 0; j <= options_.nsublevels - 1; j++) {
TEvolution aux;
aux.Lx = cv::Mat::zeros(options.img_height, options.img_width, CV_32F);
aux.Ly = cv::Mat::zeros(options.img_height, options.img_width, CV_32F);
aux.Lxx = cv::Mat::zeros(options.img_height, options.img_width, CV_32F);
aux.Lxy = cv::Mat::zeros(options.img_height, options.img_width, CV_32F);
aux.Lyy = cv::Mat::zeros(options.img_height, options.img_width, CV_32F);
aux.Lflow = cv::Mat::zeros(options.img_height, options.img_width, CV_32F);
aux.Lt = cv::Mat::zeros(options.img_height, options.img_width, CV_32F);
aux.Lsmooth = cv::Mat::zeros(options.img_height, options.img_width, CV_32F);
aux.Lstep = cv::Mat::zeros(options.img_height, options.img_width, CV_32F);
aux.Ldet = cv::Mat::zeros(options.img_height, options.img_width, CV_32F);
aux.esigma = options.soffset*pow((float)2.0f, (float)(j) / (float)(options.nsublevels)+i);
aux.Lx = cv::Mat::zeros(options_.img_height, options_.img_width, CV_32F);
aux.Ly = cv::Mat::zeros(options_.img_height, options_.img_width, CV_32F);
aux.Lxx = cv::Mat::zeros(options_.img_height, options_.img_width, CV_32F);
aux.Lxy = cv::Mat::zeros(options_.img_height, options_.img_width, CV_32F);
aux.Lyy = cv::Mat::zeros(options_.img_height, options_.img_width, CV_32F);
aux.Lt = cv::Mat::zeros(options_.img_height, options_.img_width, CV_32F);
aux.Lsmooth = cv::Mat::zeros(options_.img_height, options_.img_width, CV_32F);
aux.Ldet = cv::Mat::zeros(options_.img_height, options_.img_width, CV_32F);
aux.esigma = options_.soffset*pow((float)2.0f, (float)(j) / (float)(options_.nsublevels)+i);
aux.etime = 0.5f*(aux.esigma*aux.esigma);
aux.sigma_size = fRound(aux.esigma);
aux.octave = (float)i;
aux.sublevel = (float)j;
aux.octave = i;
aux.sublevel = j;
evolution_.push_back(aux);
}
}
// Allocate memory for the FED number of cycles and time steps
if (options.use_fed) {
for (size_t i = 1; i < evolution_.size(); i++) {
int naux = 0;
vector<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
......@@ -121,52 +98,47 @@ int KAZEFeatures::Create_Nonlinear_Scale_Space(const cv::Mat &img)
// Copy the original image to the first level of the evolution
img.copyTo(evolution_[0].Lt);
gaussian_2D_convolution(evolution_[0].Lt, evolution_[0].Lt, 0, 0, options.soffset);
gaussian_2D_convolution(evolution_[0].Lt, evolution_[0].Lsmooth, 0, 0, options.sderivatives);
gaussian_2D_convolution(evolution_[0].Lt, evolution_[0].Lt, 0, 0, options_.soffset);
gaussian_2D_convolution(evolution_[0].Lt, evolution_[0].Lsmooth, 0, 0, options_.sderivatives);
// Firstly compute the kcontrast factor
Compute_KContrast(evolution_[0].Lt, options.kcontrast_percentille);
Compute_KContrast(evolution_[0].Lt, options_.kcontrast_percentille);
// Allocate memory for the flow and step images
cv::Mat Lflow = cv::Mat::zeros(evolution_[0].Lt.rows, evolution_[0].Lt.cols, CV_32F);
cv::Mat Lstep = cv::Mat::zeros(evolution_[0].Lt.rows, evolution_[0].Lt.cols, CV_32F);
// Now generate the rest of evolution levels
for (size_t i = 1; i < evolution_.size(); i++) {
evolution_[i - 1].Lt.copyTo(evolution_[i].Lt);
gaussian_2D_convolution(evolution_[i - 1].Lt, evolution_[i].Lsmooth, 0, 0, options.sderivatives);
gaussian_2D_convolution(evolution_[i - 1].Lt, evolution_[i].Lsmooth, 0, 0, options_.sderivatives);
// Compute the Gaussian derivatives Lx and Ly
Scharr(evolution_[i].Lsmooth, evolution_[i].Lx, CV_32F, 1, 0, 1, 0, BORDER_DEFAULT);
Scharr(evolution_[i].Lsmooth, evolution_[i].Ly, CV_32F, 0, 1, 1, 0, BORDER_DEFAULT);
// Compute the conductivity equation
if (options.diffusivity == KAZEOptions::PM_G1) {
pm_g1(evolution_[i].Lx, evolution_[i].Ly, evolution_[i].Lflow, options.kcontrast);
if (options_.diffusivity == cv::DIFF_PM_G1) {
pm_g1(evolution_[i].Lx, evolution_[i].Ly, Lflow, options_.kcontrast);
}
else if (options.diffusivity == KAZEOptions::PM_G2) {
pm_g2(evolution_[i].Lx, evolution_[i].Ly, evolution_[i].Lflow, options.kcontrast);
else if (options_.diffusivity == cv::DIFF_PM_G2) {
pm_g2(evolution_[i].Lx, evolution_[i].Ly, Lflow, options_.kcontrast);
}
else if (options.diffusivity == KAZEOptions::WEICKERT) {
weickert_diffusivity(evolution_[i].Lx, evolution_[i].Ly, evolution_[i].Lflow, options.kcontrast);
else if (options_.diffusivity == cv::DIFF_WEICKERT) {
weickert_diffusivity(evolution_[i].Lx, evolution_[i].Ly, Lflow, options_.kcontrast);
}
// Perform FED n inner steps
if (options.use_fed) {
for (int j = 0; j < nsteps_[i - 1]; j++) {
nld_step_scalar(evolution_[i].Lt, evolution_[i].Lflow, evolution_[i].Lstep, tsteps_[i - 1][j]);
}
}
else {
// Perform the evolution step with AOS
AOS_Step_Scalar(evolution_[i].Lt, evolution_[i - 1].Lt, evolution_[i].Lflow,
evolution_[i].etime - evolution_[i - 1].etime);
for (int j = 0; j < nsteps_[i - 1]; j++) {
nld_step_scalar(evolution_[i].Lt, Lflow, Lstep, tsteps_[i - 1][j]);
}
}
return 0;
}
//*************************************************************************************
//*************************************************************************************
/* ************************************************************************* */
/**
* @brief This method computes the k contrast factor
* @param img Input image
......@@ -174,38 +146,10 @@ int KAZEFeatures::Create_Nonlinear_Scale_Space(const cv::Mat &img)
*/
void KAZEFeatures::Compute_KContrast(const cv::Mat &img, const float &kpercentile)
{
options.kcontrast = compute_k_percentile(img, kpercentile, options.sderivatives, options.kcontrast_bins, 0, 0);
}
//*************************************************************************************
//*************************************************************************************
/**
* @brief This method computes the multiscale derivatives for the nonlinear scale space
*/
void KAZEFeatures::Compute_Multiscale_Derivatives(void)
{
// TODO: use cv::parallel_for_
for (size_t i = 0; i < evolution_.size(); i++)
{
// Compute multiscale derivatives for the detector
compute_scharr_derivatives(evolution_[i].Lsmooth, evolution_[i].Lx, 1, 0, evolution_[i].sigma_size);
compute_scharr_derivatives(evolution_[i].Lsmooth, evolution_[i].Ly, 0, 1, evolution_[i].sigma_size);
compute_scharr_derivatives(evolution_[i].Lx, evolution_[i].Lxx, 1, 0, evolution_[i].sigma_size);
compute_scharr_derivatives(evolution_[i].Ly, evolution_[i].Lyy, 0, 1, evolution_[i].sigma_size);
compute_scharr_derivatives(evolution_[i].Lx, evolution_[i].Lxy, 0, 1, evolution_[i].sigma_size);
evolution_[i].Lx = evolution_[i].Lx*((evolution_[i].sigma_size));
evolution_[i].Ly = evolution_[i].Ly*((evolution_[i].sigma_size));
evolution_[i].Lxx = evolution_[i].Lxx*((evolution_[i].sigma_size)*(evolution_[i].sigma_size));
evolution_[i].Lxy = evolution_[i].Lxy*((evolution_[i].sigma_size)*(evolution_[i].sigma_size));
evolution_[i].Lyy = evolution_[i].Lyy*((evolution_[i].sigma_size)*(evolution_[i].sigma_size));
}
options_.kcontrast = compute_k_percentile(img, kpercentile, options_.sderivatives, options_.kcontrast_bins, 0, 0);
}
//*************************************************************************************
//*************************************************************************************
/* ************************************************************************* */
/**
* @brief This method computes the feature detector response for the nonlinear scale space
* @note We use the Hessian determinant as feature detector
......@@ -219,9 +163,9 @@ void KAZEFeatures::Compute_Detector_Response(void)
for (size_t i = 0; i < evolution_.size(); i++)
{
for (int ix = 0; ix < options.img_height; ix++)
for (int ix = 0; ix < options_.img_height; ix++)
{
for (int jx = 0; jx < options.img_width; jx++)
for (int jx = 0; jx < options_.img_width; jx++)
{
lxx = *(evolution_[i].Lxx.ptr<float>(ix)+jx);
lxy = *(evolution_[i].Lxy.ptr<float>(ix)+jx);
......@@ -232,9 +176,7 @@ void KAZEFeatures::Compute_Detector_Response(void)
}
}
//*************************************************************************************
//*************************************************************************************
/* ************************************************************************* */
/**
* @brief This method selects interesting keypoints through the nonlinear scale space
* @param kpts Vector of keypoints
......@@ -242,27 +184,127 @@ void KAZEFeatures::Compute_Detector_Response(void)
void KAZEFeatures::Feature_Detection(std::vector<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;
......@@ -283,10 +325,8 @@ void KAZEFeatures::Determinant_Hessian_Parallel(std::vector<cv::KeyPoint>& kpts)
kpts_par_.push_back(aux);
}
// TODO: Use cv::parallel_for_
for (int i = 1; i < (int)evolution_.size() - 1; i++) {
Find_Extremum_Threading(i);
}
cv::parallel_for_(cv::Range(1, (int)evolution_.size()-1),
FindExtremumKAZEInvoker(evolution_, kpts_par_, options_));
// Now fill the vector of keypoints!!!
for (int i = 0; i < (int)kpts_par_.size(); i++) {
......@@ -343,63 +383,7 @@ void KAZEFeatures::Determinant_Hessian_Parallel(std::vector<cv::KeyPoint>& kpts)
}
}
//*************************************************************************************
//*************************************************************************************
/**
* @brief This method is called by the thread which is responsible of finding extrema
* at a given nonlinear scale level
* @param level Index in the nonlinear scale space evolution
*/
void KAZEFeatures::Find_Extremum_Threading(const int& level) {
float value = 0.0;
bool is_extremum = false;
for (int ix = 1; ix < options.img_height - 1; ix++) {
for (int jx = 1; jx < options.img_width - 1; jx++) {
is_extremum = false;
value = *(evolution_[level].Ldet.ptr<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
......@@ -475,10 +459,10 @@ void KAZEFeatures::Do_Subpixel_Refinement(std::vector<cv::KeyPoint> &kpts) {
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
......@@ -497,17 +481,15 @@ void KAZEFeatures::Do_Subpixel_Refinement(std::vector<cv::KeyPoint> &kpts) {
}
}
//*************************************************************************************
//*************************************************************************************
/* ************************************************************************* */
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)
{
}
......@@ -517,26 +499,26 @@ public:
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));
......@@ -549,12 +531,13 @@ private:
void Get_KAZE_Upright_Descriptor_128(const cv::KeyPoint& kpt, float* desc) const;
void Get_KAZE_Descriptor_128(const cv::KeyPoint& kpt, float *desc) const;
std::vector<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
......@@ -562,20 +545,23 @@ private:
*/
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
......@@ -651,9 +637,7 @@ void KAZEFeatures::Compute_Main_Orientation(cv::KeyPoint &kpt, const std::vector
}
}
//*************************************************************************************
//*************************************************************************************
/* ************************************************************************* */
/**
* @brief This method computes the upright descriptor (not rotation invariant) of
* the provided keypoint
......@@ -673,7 +657,7 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Upright_Descriptor_64(const cv::KeyPoint
float fx = 0.0, fy = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0, res4 = 0.0;
int dsize = 0, scale = 0, level = 0;
std::vector<TEvolution>& evolution_ = *_evolution;
std::vector<TEvolution>& evolution = *evolution_;
// Subregion centers for the 4x4 gaussian weighting
float cx = -0.5f, cy = 0.5f;
......@@ -724,26 +708,26 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Upright_Descriptor_64(const cv::KeyPoint
y1 = (int)(sample_y - 0.5f);
x1 = (int)(sample_x - 0.5f);
checkDescriptorLimits(x1, y1, options.img_width, options.img_height);
checkDescriptorLimits(x1, y1, options_.img_width, options_.img_height);
y2 = (int)(sample_y + 0.5f);
x2 = (int)(sample_x + 0.5f);
checkDescriptorLimits(x2, y2, options.img_width, options.img_height);
checkDescriptorLimits(x2, y2, options_.img_width, options_.img_height);
fx = sample_x - x1;
fy = sample_y - y1;
res1 = *(evolution_[level].Lx.ptr<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;
......@@ -779,15 +763,9 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Upright_Descriptor_64(const cv::KeyPoint
for (i = 0; i < dsize; i++) {
desc[i] /= len;
}
if (options.use_clipping_normalilzation) {
clippingDescriptor(desc, dsize, options.clipping_normalization_niter, options.clipping_normalization_ratio);
}
}
//*************************************************************************************
//*************************************************************************************
/* ************************************************************************* */
/**
* @brief This method computes the descriptor of the provided keypoint given the
* main orientation of the keypoint
......@@ -807,7 +785,7 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Descriptor_64(const cv::KeyPoint &kpt, fl
int kx = 0, ky = 0, i = 0, j = 0, dcount = 0;
int dsize = 0, scale = 0, level = 0;
std::vector<TEvolution>& evolution_ = *_evolution;
std::vector<TEvolution>& evolution = *evolution_;
// Subregion centers for the 4x4 gaussian weighting
float cx = -0.5f, cy = 0.5f;
......@@ -862,26 +840,26 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Descriptor_64(const cv::KeyPoint &kpt, fl
y1 = fRound(sample_y - 0.5f);
x1 = fRound(sample_x - 0.5f);
checkDescriptorLimits(x1, y1, options.img_width, options.img_height);
checkDescriptorLimits(x1, y1, options_.img_width, options_.img_height);
y2 = (int)(sample_y + 0.5f);
x2 = (int)(sample_x + 0.5f);
checkDescriptorLimits(x2, y2, options.img_width, options.img_height);
checkDescriptorLimits(x2, y2, options_.img_width, options_.img_height);
fx = sample_x - x1;
fy = sample_y - y1;
res1 = *(evolution_[level].Lx.ptr<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
......@@ -914,15 +892,9 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Descriptor_64(const cv::KeyPoint &kpt, fl
for (i = 0; i < dsize; i++) {
desc[i] /= len;
}
if (options.use_clipping_normalilzation) {
clippingDescriptor(desc, dsize, options.clipping_normalization_niter, options.clipping_normalization_ratio);
}
}
//*************************************************************************************
//*************************************************************************************
/* ************************************************************************* */
/**
* @brief This method computes the extended upright descriptor (not rotation invariant) of
* the provided keypoint
......@@ -947,7 +919,7 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Upright_Descriptor_128(const cv::KeyPoint
// Subregion centers for the 4x4 gaussian weighting
float cx = -0.5f, cy = 0.5f;
std::vector<TEvolution>& evolution_ = *_evolution;
std::vector<TEvolution>& evolution = *evolution_;
// Set the descriptor size and the sample and pattern sizes
dsize = 128;
......@@ -998,26 +970,26 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Upright_Descriptor_128(const cv::KeyPoint
y1 = (int)(sample_y - 0.5f);
x1 = (int)(sample_x - 0.5f);
checkDescriptorLimits(x1, y1, options.img_width, options.img_height);
checkDescriptorLimits(x1, y1, options_.img_width, options_.img_height);
y2 = (int)(sample_y + 0.5f);
x2 = (int)(sample_x + 0.5f);
checkDescriptorLimits(x2, y2, options.img_width, options.img_height);
checkDescriptorLimits(x2, y2, options_.img_width, options_.img_height);
fx = sample_x - x1;
fy = sample_y - y1;
res1 = *(evolution_[level].Lx.ptr<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;
......@@ -1072,15 +1044,9 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Upright_Descriptor_128(const cv::KeyPoint
for (i = 0; i < dsize; i++) {
desc[i] /= len;
}
if (options.use_clipping_normalilzation) {
clippingDescriptor(desc, dsize, options.clipping_normalization_niter, options.clipping_normalization_ratio);
}
}
//*************************************************************************************
//*************************************************************************************
/* ************************************************************************* */
/**
* @brief This method computes the extended G-SURF descriptor of the provided keypoint
* given the main orientation of the keypoint
......@@ -1102,7 +1068,7 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Descriptor_128(const cv::KeyPoint &kpt, f
int kx = 0, ky = 0, i = 0, j = 0, dcount = 0;
int dsize = 0, scale = 0, level = 0;
std::vector<TEvolution>& evolution_ = *_evolution;
std::vector<TEvolution>& evolution = *evolution_;
// Subregion centers for the 4x4 gaussian weighting
float cx = -0.5f, cy = 0.5f;
......@@ -1160,26 +1126,26 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Descriptor_128(const cv::KeyPoint &kpt, f
y1 = fRound(sample_y - 0.5f);
x1 = fRound(sample_x - 0.5f);
checkDescriptorLimits(x1, y1, options.img_width, options.img_height);
checkDescriptorLimits(x1, y1, options_.img_width, options_.img_height);
y2 = (int)(sample_y + 0.5f);
x2 = (int)(sample_x + 0.5f);
checkDescriptorLimits(x2, y2, options.img_width, options.img_height);
checkDescriptorLimits(x2, y2, options_.img_width, options_.img_height);
fx = sample_x - x1;
fy = sample_y - y1;
res1 = *(evolution_[level].Lx.ptr<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
......@@ -1235,298 +1201,4 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Descriptor_128(const cv::KeyPoint &kpt, f
for (i = 0; i < dsize; i++) {
desc[i] /= len;
}
if (options.use_clipping_normalilzation) {
clippingDescriptor(desc, dsize, options.clipping_normalization_niter, options.clipping_normalization_ratio);
}
}
//*************************************************************************************
//*************************************************************************************
/**
* @brief This method performs a scalar non-linear diffusion step using AOS schemes
* @param Ld Image at a given evolution step
* @param Ldprev Image at a previous evolution step
* @param c Conductivity image
* @param stepsize Stepsize for the nonlinear diffusion evolution
* @note If c is constant, the diffusion will be linear
* If c is a matrix of the same size as Ld, the diffusion will be nonlinear
* The stepsize can be arbitrarily large
*/
void KAZEFeatures::AOS_Step_Scalar(cv::Mat &Ld, const cv::Mat &Ldprev, const cv::Mat &c, const float& stepsize) {
AOS_Rows(Ldprev, c, stepsize);
AOS_Columns(Ldprev, c, stepsize);
Ld = 0.5f*(Lty_ + Ltx_.t());
}
//*************************************************************************************
//*************************************************************************************
/**
* @brief This method performs performs 1D-AOS for the image rows
* @param Ldprev Image at a previous evolution step
* @param c Conductivity image
* @param stepsize Stepsize for the nonlinear diffusion evolution
*/
void KAZEFeatures::AOS_Rows(const cv::Mat &Ldprev, const cv::Mat &c, const float& stepsize) {
// Operate on rows
for (int i = 0; i < qr_.rows; i++) {
for (int j = 0; j < qr_.cols; j++) {
*(qr_.ptr<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);
}
modules/features2d/src/kaze/KAZEFeatures.h
View file @ 7f829608
......@@ -7,84 +7,53 @@
* @author Pablo F. Alcantarilla
*/
#ifndef KAZE_H_
#define KAZE_H_
//*************************************************************************************
//*************************************************************************************
#ifndef __OPENCV_FEATURES_2D_KAZE_FEATURES_H__
#define __OPENCV_FEATURES_2D_KAZE_FEATURES_H__
/* ************************************************************************* */
// Includes
#include "KAZEConfig.h"
#include "nldiffusion_functions.h"
#include "fed.h"
#include "TEvolution.h"
//*************************************************************************************
//*************************************************************************************
/* ************************************************************************* */
// KAZE Class Declaration
class KAZEFeatures {
private:
KAZEOptions options;
// Parameters of the Nonlinear diffusion class
std::vector<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
modules/features2d/src/kaze/TEvolution.h 0 → 100644
View file @ 7f829608
/**
* @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
modules/features2d/src/kaze/fed.h
View file @ 7f829608
#ifndef FED_H
#define FED_H
#ifndef __OPENCV_FEATURES_2D_FED_H__
#define __OPENCV_FEATURES_2D_FED_H__
//******************************************************************************
//******************************************************************************
......@@ -22,4 +22,4 @@ bool fed_is_prime_internal(const int& number);
//*************************************************************************************
//*************************************************************************************
#endif // FED_H
#endif // __OPENCV_FEATURES_2D_FED_H__
modules/features2d/src/kaze/nldiffusion_functions.h
View file @ 7f829608
......@@ -8,8 +8,8 @@
* @author Pablo F. Alcantarilla
*/
#ifndef KAZE_NLDIFFUSION_FUNCTIONS_H
#define KAZE_NLDIFFUSION_FUNCTIONS_H
#ifndef __OPENCV_FEATURES_2D_NLDIFFUSION_FUNCTIONS_H__
#define __OPENCV_FEATURES_2D_NLDIFFUSION_FUNCTIONS_H__
/* ************************************************************************* */
// Includes
......
modules/features2d/src/kaze/utils.h 0 → 100644
View file @ 7f829608
#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
modules/features2d/test/test_descriptors_regression.cpp
View file @ 7f829608
......@@ -101,8 +101,14 @@ public:
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() {}
......@@ -189,7 +195,6 @@ protected:
// Read the test image.
string imgFilename = string(ts->get_data_path()) + FEATURES2D_DIR + "/" + IMAGE_FILENAME;
Mat img = imread( imgFilename );
if( img.empty() )
{
......@@ -197,13 +202,15 @@ protected:
ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_TEST_DATA );
return;
}
vector<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 );
......@@ -244,7 +251,7 @@ protected:
}
}
}
else
if(!fs.isOpened())
{
ts->printf( cvtest::TS::LOG, "Compute and write keypoints.\n" );
fs.open( string(ts->get_data_path()) + FEATURES2D_DIR + "/keypoints.xml.gz", FileStorage::WRITE );
......@@ -295,6 +302,7 @@ protected:
const DistanceType maxDist;
Ptr<DescriptorExtractor> dextractor;
Distance distance;
Ptr<FeatureDetector> detector;
private:
CV_DescriptorExtractorTest& operator=(const CV_DescriptorExtractorTest&) { return *this; }
......@@ -340,3 +348,19 @@ TEST( Features2d_DescriptorExtractor_OpponentBRIEF, regression )
DescriptorExtractor::create("OpponentBRIEF") );
test.safe_run();
}
TEST( Features2d_DescriptorExtractor_KAZE, regression )
{
CV_DescriptorExtractorTest< L2<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();
}
modules/features2d/test/test_detectors_regression.cpp
View file @ 7f829608
......@@ -289,6 +289,18 @@ TEST( Features2d_Detector_ORB, regression )
test.safe_run();
}
TEST( Features2d_Detector_KAZE, regression )
{
CV_FeatureDetectorTest test( "detector-kaze", FeatureDetector::create("KAZE") );
test.safe_run();
}
TEST( Features2d_Detector_AKAZE, regression )
{
CV_FeatureDetectorTest test( "detector-akaze", FeatureDetector::create("AKAZE") );
test.safe_run();
}
TEST( Features2d_Detector_GridFAST, regression )
{
CV_FeatureDetectorTest test( "detector-grid-fast", FeatureDetector::create("GridFAST") );
......
modules/features2d/test/test_keypoints.cpp
View file @ 7f829608
......@@ -167,19 +167,17 @@ TEST(Features2d_Detector_Keypoints_Dense, validation)
test.safe_run();
}
// FIXIT #2807 Crash on Windows 7 x64 MSVS 2012, Linux Fedora 19 x64 with GCC 4.8.2, Linux Ubuntu 14.04 LTS x64 with GCC 4.8.2
TEST(Features2d_Detector_Keypoints_KAZE, DISABLED_validation)
TEST(Features2d_Detector_Keypoints_KAZE, validation)
{
CV_FeatureDetectorKeypointsTest test(Algorithm::create<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();
}
modules/features2d/test/test_rotation_and_scale_invariance.cpp
View file @ 7f829608
......@@ -652,8 +652,7 @@ TEST(Features2d_ScaleInvariance_Detector_BRISK, regression)
test.safe_run();
}
// FIXIT #2807 Crash on Windows 7 x64 MSVS 2012, Linux Fedora 19 x64 with GCC 4.8.2, Linux Ubuntu 14.04 LTS x64 with GCC 4.8.2
TEST(Features2d_ScaleInvariance_Detector_KAZE, DISABLED_regression)
TEST(Features2d_ScaleInvariance_Detector_KAZE, regression)
{
DetectorScaleInvarianceTest test(Algorithm::create<FeatureDetector>("Feature2D.KAZE"),
0.08f,
......@@ -661,8 +660,7 @@ TEST(Features2d_ScaleInvariance_Detector_KAZE, DISABLED_regression)
test.safe_run();
}
// FIXIT #2807 Crash on Windows 7 x64 MSVS 2012, Linux Fedora 19 x64 with GCC 4.8.2, Linux Ubuntu 14.04 LTS x64 with GCC 4.8.2
TEST(Features2d_ScaleInvariance_Detector_AKAZE, DISABLED_regression)
TEST(Features2d_ScaleInvariance_Detector_AKAZE, regression)
{
DetectorScaleInvarianceTest test(Algorithm::create<FeatureDetector>("Feature2D.AKAZE"),
0.08f,
......
samples/cpp/H1to3p.xml 0 → 100755
View file @ 7f829608
<?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>
samples/cpp/tutorial_code/features2D/AKAZE_match.cpp 0 → 100755
View file @ 7f829608
#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;
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment