KAZE and AKAZE integration initial commit

modules/features2d/src/akaze/AKAZE.h

@@ -0,0 +1,175 @@
+/**
+ * @file AKAZE.h
+ * @brief Main class for detecting and computing binary descriptors in an
+ * accelerated nonlinear scale space
+ * @date Mar 27, 2013
+ * @author Pablo F. Alcantarilla, Jesus Nuevo
+ */
+
+#ifndef _AKAZE_H_
+#define _AKAZE_H_
+
+//*************************************************************************************
+//*************************************************************************************
+
+// Includes
+#include "config.h"
+#include "fed.h"
+#include "utils.h"
+#include "nldiffusion_functions.h"
+
+//*************************************************************************************
+//*************************************************************************************
+
+// AKAZE Class Declaration
+class AKAZE {
+
+private:
+
+  // Parameters of the AKAZE class
+  int omax_;               // Maximum octave level
+  int noctaves_;           // Number of octaves
+  int nsublevels_;         // Number of sublevels per octave level
+  int img_width_;          // Width of the original image
+  int img_height_;         // Height of the original image
+  float soffset_;          // Base scale offset
+  float factor_size_;      // Factor for the multiscale derivatives
+  float sderivatives_;     // Standard deviation of the Gaussian for the nonlinear diff. derivatives
+  float kcontrast_;        // The contrast parameter for the scalar nonlinear diffusion
+  float dthreshold_;       // Feature detector threshold response
+  int diffusivity_;        // Diffusivity type, 0->PM G1, 1->PM G2, 2-> Weickert, 3->Charbonnier
+  int descriptor_;         // Descriptor mode:
+                           // 0-> SURF_UPRIGHT, 1->SURF
+                           // 2-> M-SURF_UPRIGHT, 3->M-SURF
+                           // 4-> M-LDB_UPRIGHT, 5->M-LDB
+  int descriptor_size_;    // Size of the descriptor in bits. Use 0 for the full descriptor
+  int descriptor_pattern_size_;    // Size of the pattern. Actual size sampled is 2*pattern_size
+  int descriptor_channels_;        // Number of channels to consider in the M-LDB descriptor
+  bool save_scale_space_; // For saving scale space images
+  bool verbosity_;	// Verbosity level
+  std::vector<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
+
+  // Some matrices for the M-LDB descriptor computation
+  cv::Mat descriptorSamples_;  // List of positions in the grids to sample LDB bits from.
+  cv::Mat descriptorBits_;
+  cv::Mat bitMask_;
+
+  // Computation times variables in ms
+  double tkcontrast_;      // Kcontrast factor computation
+  double tscale_;          // Nonlinear Scale space generation
+  double tderivatives_;    // Multiscale derivatives
+  double tdetector_;       // Feature detector
+  double textrema_;        // Scale Space extrema
+  double tsubpixel_;       // Subpixel refinement
+  double tdescriptor_;     // Feature descriptors
+
+public:
+
+  // Constructor
+  AKAZE(const AKAZEOptions &options);
+
+  // Destructor
+  ~AKAZE(void);
+
+  // Setters
+  void Set_Octave_Max(const int& omax) {
+    omax_ = omax;
+  }
+  void Set_NSublevels(const int& nsublevels) {
+    nsublevels_ = nsublevels;
+  }
+  void Set_Save_Scale_Space_Flag(const bool& save_scale_space) {
+    save_scale_space_ = save_scale_space;
+  }
+  void Set_Image_Width(const int& img_width) {
+    img_width_ = img_width;
+  }
+  void Set_Image_Height(const int& img_height) {
+    img_height_ = img_height;
+  }
+
+  // Getters
+  int Get_Image_Width(void) {
+    return img_width_;
+  }
+  int Get_Image_Height(void) {
+    return img_height_;
+  }
+  double Get_Time_KContrast(void) {
+    return tkcontrast_;
+  }
+  double Get_Time_Scale_Space(void) {
+    return tscale_;
+  }
+  double Get_Time_Derivatives(void) {
+    return tderivatives_;
+  }
+  double Get_Time_Detector(void) {
+    return tdetector_;
+  }
+  double Get_Time_Descriptor(void) {
+    return tdescriptor_;
+  }
+
+  // Scale Space methods
+  void Allocate_Memory_Evolution(void);
+  int Create_Nonlinear_Scale_Space(const cv::Mat& img);
+  void Feature_Detection(std::vector<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);
+  void Feature_Suppression_Distance(std::vector<cv::KeyPoint>& kpts, float mdist);
+
+  // Feature description methods
+  void Compute_Descriptors(std::vector<cv::KeyPoint>& kpts, cv::Mat& desc);
+  void Compute_Main_Orientation_SURF(cv::KeyPoint& kpt);
+
+  // SURF Pattern Descriptor
+  void Get_SURF_Descriptor_Upright_64(const cv::KeyPoint& kpt, float *desc);
+  void Get_SURF_Descriptor_64(const cv::KeyPoint& kpt, float *desc);
+
+  // M-SURF Pattern Descriptor
+  void Get_MSURF_Upright_Descriptor_64(const cv::KeyPoint& kpt, float *desc);
+  void Get_MSURF_Descriptor_64(const cv::KeyPoint& kpt, float *desc);
+
+  // M-LDB Pattern Descriptor
+  void Get_Upright_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned char *desc);
+  void Get_MLDB_Full_Descriptor(const cv::KeyPoint& kpt, unsigned char *desc);
+  void Get_Upright_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char *desc);
+  void Get_MLDB_Descriptor_Subset(const cv::KeyPoint& kpt, unsigned char *desc);
+
+  // Methods for saving some results and showing computation times
+  void Save_Scale_Space(void);
+  void Save_Detector_Responses(void);
+  void Show_Computation_Times(void);
+};
+
+//*************************************************************************************
+//*************************************************************************************
+
+// Inline functions
+/**
+ * @brief This function sets default parameters for the A-KAZE detector.
+ * @param options AKAZE options
+ */
+void setDefaultAKAZEOptions(AKAZEOptions& options);
+
+// Inline functions
+void generateDescriptorSubsample(cv::Mat& sampleList, cv::Mat& comparisons,
+                                 int nbits, int pattern_size, int nchannels);
+float get_angle(float x, float y);
+float gaussian(float x, float y, float sigma);
+void check_descriptor_limits(int& x, int& y, const int width, const int height);
+int fRound(float flt);
+
+//*************************************************************************************
+//*************************************************************************************
+
+#endif

modules/features2d/src/akaze/config.h

@@ -0,0 +1,155 @@
+#ifndef _CONFIG_H_
+#define _CONFIG_H_
+
+// STL
+#include <string>
+#include <vector>
+#include <cmath>
+#include <bitset>
+#include <iomanip>
+
+// OpenCV
+#include "precomp.hpp"
+
+// OpenMP
+#ifdef _OPENMP
+# include <omp.h>
+#endif
+
+// Lookup table for 2d gaussian (sigma = 2.5) where (0,0) is top left and (6,6) is bottom right
+const float gauss25[7][7] = {
+  {0.02546481f,	0.02350698f,	0.01849125f,	0.01239505f,	0.00708017f,	0.00344629f,	0.00142946f},
+  {0.02350698f,	0.02169968f,	0.01706957f,	0.01144208f,	0.00653582f,	0.00318132f,	0.00131956f},
+  {0.01849125f,	0.01706957f,	0.01342740f,	0.00900066f,	0.00514126f,	0.00250252f,	0.00103800f},
+  {0.01239505f,	0.01144208f,	0.00900066f,	0.00603332f,	0.00344629f,	0.00167749f,	0.00069579f},
+  {0.00708017f,	0.00653582f,	0.00514126f,	0.00344629f,	0.00196855f,	0.00095820f,	0.00039744f},
+  {0.00344629f,	0.00318132f,	0.00250252f,	0.00167749f,	0.00095820f,	0.00046640f,	0.00019346f},
+  {0.00142946f,	0.00131956f,	0.00103800f,	0.00069579f,	0.00039744f,	0.00019346f,	0.00008024f}
+};
+
+
+// Scale Space parameters
+const float DEFAULT_SCALE_OFFSET = 1.60f;    // Base scale offset (sigma units)
+const float DEFAULT_FACTOR_SIZE = 1.5f;      // Factor for the multiscale derivatives
+const int DEFAULT_OCTAVE_MIN = 0;            // Initial octave level (-1 means that the size of the input image is duplicated)
+const int DEFAULT_OCTAVE_MAX = 4;            // Maximum octave evolution of the image 2^sigma (coarsest scale sigma units)
+const int DEFAULT_NSUBLEVELS = 4;            // Default number of sublevels per scale level
+const int DEFAULT_DIFFUSIVITY_TYPE = 1;
+const float KCONTRAST_PERCENTILE = 0.7f;
+const int KCONTRAST_NBINS = 300;
+const float DEFAULT_SIGMA_SMOOTHING_DERIVATIVES = 1.0f;
+const float DEFAULT_KCONTRAST = .01f;
+
+
+// Detector Parameters
+const float DEFAULT_DETECTOR_THRESHOLD = 0.001f;           // Detector response threshold to accept point
+const float DEFAULT_MIN_DETECTOR_THRESHOLD = 0.00001f;     // Minimum Detector response threshold to accept point
+const int DEFAULT_LDB_DESCRIPTOR_SIZE = 0;  // Use 0 for the full descriptor, or the number of bits
+const int DEFAULT_LDB_PATTERN_SIZE = 10;    // Actual patch size is 2*pattern_size*point.scale;
+const int DEFAULT_LDB_CHANNELS = 3;
+
+// Descriptor Parameters
+enum DESCRIPTOR_TYPE
+{
+  SURF_UPRIGHT = 0, // Upright descriptors, not invariant to rotation
+  SURF = 1,
+  MSURF_UPRIGHT = 2, // Upright descriptors, not invariant to rotation
+  MSURF = 3,
+  MLDB_UPRIGHT = 4, // Upright descriptors, not invariant to rotation
+  MLDB = 5
+};
+
+const int DEFAULT_DESCRIPTOR = MLDB;
+
+// Some debugging options
+const bool DEFAULT_SAVE_SCALE_SPACE = false; // For saving the scale space images
+const bool DEFAULT_VERBOSITY = false; // Verbosity level (0->no verbosity)
+const bool DEFAULT_SHOW_RESULTS = true; // For showing the output image with the detected features plus some ratios
+const bool DEFAULT_SAVE_KEYPOINTS = false; // For saving the list of keypoints
+
+// Options structure
+struct AKAZEOptions
+{
+  int omin;
+  int omax;
+  int nsublevels;
+  int img_width;
+  int img_height;
+  int diffusivity;
+  float soffset;
+  float sderivatives;
+  float dthreshold;
+  float dthreshold2;
+  int descriptor;
+  int descriptor_size;
+  int descriptor_channels;
+  int descriptor_pattern_size;
+  bool save_scale_space;
+  bool save_keypoints;
+  bool verbosity;
+
+  AKAZEOptions()
+  {
+    // Load the default options
+    soffset = DEFAULT_SCALE_OFFSET;
+    omax = DEFAULT_OCTAVE_MAX;
+    nsublevels = DEFAULT_NSUBLEVELS;
+    dthreshold = DEFAULT_DETECTOR_THRESHOLD;
+    diffusivity = DEFAULT_DIFFUSIVITY_TYPE;
+    descriptor = DEFAULT_DESCRIPTOR;
+    descriptor_size = DEFAULT_LDB_DESCRIPTOR_SIZE;
+    descriptor_channels = DEFAULT_LDB_CHANNELS;
+    descriptor_pattern_size = DEFAULT_LDB_PATTERN_SIZE;
+    sderivatives = DEFAULT_SIGMA_SMOOTHING_DERIVATIVES;
+    save_scale_space = DEFAULT_SAVE_SCALE_SPACE;
+    save_keypoints = DEFAULT_SAVE_KEYPOINTS;
+    verbosity = DEFAULT_VERBOSITY;
+  }
+
+  friend std::ostream& operator<<(std::ostream& os,
+                                  const AKAZEOptions& akaze_options)
+  {
+    os << std::left;
+#define CHECK_AKAZE_OPTION(option) \
+    os << std::setw(33) << #option << " =  " << option << std::endl
+
+    // Scale-space parameters.
+    CHECK_AKAZE_OPTION(akaze_options.omax);
+    CHECK_AKAZE_OPTION(akaze_options.nsublevels);
+    CHECK_AKAZE_OPTION(akaze_options.soffset);
+    CHECK_AKAZE_OPTION(akaze_options.sderivatives);
+    CHECK_AKAZE_OPTION(akaze_options.diffusivity);
+    // Detection parameters.
+    CHECK_AKAZE_OPTION(akaze_options.dthreshold);
+    // Descriptor parameters.
+    CHECK_AKAZE_OPTION(akaze_options.descriptor);
+    CHECK_AKAZE_OPTION(akaze_options.descriptor_channels);
+    CHECK_AKAZE_OPTION(akaze_options.descriptor_size);
+    // Save scale-space
+    CHECK_AKAZE_OPTION(akaze_options.save_scale_space);
+    // Verbose option for debug.
+    CHECK_AKAZE_OPTION(akaze_options.verbosity);
+#undef CHECK_AKAZE_OPTIONS
+
+    return os;
+  }
+};
+
+struct tevolution
+{
+	cv::Mat Lx, Ly;	// First order spatial derivatives
+	cv::Mat Lxx, Lxy, Lyy;	// Second order spatial derivatives
+	cv::Mat Lflow;	// Diffusivity image
+	cv::Mat Lt;	// Evolution image
+	cv::Mat Lsmooth; // Smoothed image
+	cv::Mat Lstep; // Evolution step update
+	cv::Mat Ldet; // Detector response
+	float etime;	// Evolution time
+	float esigma;	// Evolution sigma. For linear diffusion t = sigma^2 / 2
+  int octave;	// Image octave
+  int sublevel;	// Image sublevel in each octave
+	int sigma_size;	// Integer sigma. For computing the feature detector responses
+};
+
+
+#endif

modules/features2d/src/akaze/fed.h

@@ -0,0 +1,26 @@
+#ifndef FED_H
+#define FED_H
+
+//******************************************************************************
+//******************************************************************************
+
+// Includes
+#include <iostream>
+#include <vector>
+
+//*************************************************************************************
+//*************************************************************************************
+
+// Declaration of functions
+int fed_tau_by_process_time(const float& T, const int& M, const float& tau_max,
+                            const bool& reordering, std::vector<float>& tau);
+int fed_tau_by_cycle_time(const float& t, const float& tau_max,
+                          const bool& reordering, std::vector<float> &tau) ;
+int fed_tau_internal(const int& n, const float& scale, const float& tau_max,
+                     const bool& reordering, std::vector<float> &tau);
+bool fed_is_prime_internal(const int& number);
+
+//*************************************************************************************
+//*************************************************************************************
+
+#endif // FED_H

modules/features2d/src/akaze/nldiffusion_functions.cpp

@@ -0,0 +1,431 @@
+//=============================================================================
+//
+// nldiffusion_functions.cpp
+// Authors: Pablo F. Alcantarilla (1), Jesus Nuevo (2)
+// Institutions: Georgia Institute of Technology (1)
+//               TrueVision Solutions (2)
+// Date: 15/09/2013
+// Email: pablofdezalc@gmail.com
+//
+// AKAZE Features Copyright 2013, Pablo F. Alcantarilla, Jesus Nuevo
+// All Rights Reserved
+// See LICENSE for the license information
+//=============================================================================
+
+/**
+ * @file nldiffusion_functions.cpp
+ * @brief Functions for nonlinear diffusion filtering applications
+ * @date Sep 15, 2013
+ * @author Pablo F. Alcantarilla, Jesus Nuevo
+ */
+
+#include "nldiffusion_functions.h"
+
+using namespace std;
+using namespace cv;
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function smoothes an image with a Gaussian kernel
+ * @param src Input image
+ * @param dst Output image
+ * @param ksize_x Kernel size in X-direction (horizontal)
+ * @param ksize_y Kernel size in Y-direction (vertical)
+ * @param sigma Kernel standard deviation
+ */
+void gaussian_2D_convolution(const cv::Mat& src, cv::Mat& dst, const size_t& ksize_x,
+                             const size_t& ksize_y, const float& sigma) {
+
+  size_t ksize_x_ = 0, ksize_y_ = 0;
+
+  // Compute an appropriate kernel size according to the specified sigma
+  if (sigma > ksize_x || sigma > ksize_y || ksize_x == 0 || ksize_y == 0) {
+    ksize_x_ = ceil(2.0*(1.0 + (sigma-0.8)/(0.3)));
+    ksize_y_ = ksize_x_;
+  }
+
+  // The kernel size must be and odd number
+  if ((ksize_x_ % 2) == 0) {
+    ksize_x_ += 1;
+  }
+
+  if ((ksize_y_ % 2) == 0) {
+    ksize_y_ += 1;
+  }
+
+  // Perform the Gaussian Smoothing with border replication
+  GaussianBlur(src,dst,Size(ksize_x_,ksize_y_),sigma,sigma,BORDER_REPLICATE);
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function computes image derivatives with Scharr kernel
+ * @param src Input image
+ * @param dst Output image
+ * @param xorder Derivative order in X-direction (horizontal)
+ * @param yorder Derivative order in Y-direction (vertical)
+ * @note Scharr operator approximates better rotation invariance than
+ * other stencils such as Sobel. See Weickert and Scharr,
+ * A Scheme for Coherence-Enhancing Diffusion Filtering with Optimized Rotation Invariance,
+ * Journal of Visual Communication and Image Representation 2002
+ */
+void image_derivatives_scharr(const cv::Mat& src, cv::Mat& dst,
+                              const size_t& xorder, const size_t& yorder) {
+  Scharr(src,dst,CV_32F,xorder,yorder,1.0,0,BORDER_DEFAULT);
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function computes the Perona and Malik conductivity coefficient g1
+ * g1 = exp(-|dL|^2/k^2)
+ * @param Lx First order image derivative in X-direction (horizontal)
+ * @param Ly First order image derivative in Y-direction (vertical)
+ * @param dst Output image
+ * @param k Contrast factor parameter
+ */
+void pm_g1(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k) {
+  exp(-(Lx.mul(Lx)+Ly.mul(Ly))/(k*k),dst);
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function computes the Perona and Malik conductivity coefficient g2
+ * g2 = 1 / (1 + dL^2 / k^2)
+ * @param Lx First order image derivative in X-direction (horizontal)
+ * @param Ly First order image derivative in Y-direction (vertical)
+ * @param dst Output image
+ * @param k Contrast factor parameter
+ */
+void pm_g2(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k) {
+  dst = 1.0/(1.0+(Lx.mul(Lx)+Ly.mul(Ly))/(k*k));
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function computes Weickert conductivity coefficient gw
+ * @param Lx First order image derivative in X-direction (horizontal)
+ * @param Ly First order image derivative in Y-direction (vertical)
+ * @param dst Output image
+ * @param k Contrast factor parameter
+ * @note For more information check the following paper: J. Weickert
+ * Applications of nonlinear diffusion in image processing and computer vision,
+ * Proceedings of Algorithmy 2000
+ */
+void weickert_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k) {
+  Mat modg;
+  pow((Lx.mul(Lx) + Ly.mul(Ly))/(k*k),4,modg);
+  cv::exp(-3.315/modg, dst);
+  dst = 1.0 - dst;
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function computes Charbonnier conductivity coefficient gc
+ * gc = 1 / sqrt(1 + dL^2 / k^2)
+ * @param Lx First order image derivative in X-direction (horizontal)
+ * @param Ly First order image derivative in Y-direction (vertical)
+ * @param dst Output image
+ * @param k Contrast factor parameter
+ * @note For more information check the following paper: J. Weickert
+ * Applications of nonlinear diffusion in image processing and computer vision,
+ * Proceedings of Algorithmy 2000
+ */
+void charbonnier_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k) {
+  Mat den;
+  cv::sqrt(1.0+(Lx.mul(Lx)+Ly.mul(Ly))/(k*k),den);
+  dst = 1.0/ den;
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function computes a good empirical value for the k contrast factor
+ * given an input image, the percentile (0-1), the gradient scale and the number of
+ * bins in the histogram
+ * @param img Input image
+ * @param perc Percentile of the image gradient histogram (0-1)
+ * @param gscale Scale for computing the image gradient histogram
+ * @param nbins Number of histogram bins
+ * @param ksize_x Kernel size in X-direction (horizontal) for the Gaussian smoothing kernel
+ * @param ksize_y Kernel size in Y-direction (vertical) for the Gaussian smoothing kernel
+ * @return k contrast factor
+ */
+float compute_k_percentile(const cv::Mat& img, const float& perc, const float& gscale,
+                           const size_t& nbins, const size_t& ksize_x, const size_t& ksize_y) {
+
+  size_t nbin = 0, nelements = 0, nthreshold = 0, k = 0;
+  float kperc = 0.0, modg = 0.0, lx = 0.0, ly = 0.0;
+  float npoints = 0.0;
+  float hmax = 0.0;
+
+  // Create the array for the histogram
+  float *hist = new float[nbins];
+
+  // Create the matrices
+  Mat gaussian = Mat::zeros(img.rows,img.cols,CV_32F);
+  Mat Lx = Mat::zeros(img.rows,img.cols,CV_32F);
+  Mat Ly = Mat::zeros(img.rows,img.cols,CV_32F);
+
+  // Set the histogram to zero, just in case
+  for (size_t i = 0; i < nbins; i++) {
+    hist[i] = 0.0;
+  }
+
+  // Perform the Gaussian convolution
+  gaussian_2D_convolution(img,gaussian,ksize_x,ksize_y,gscale);
+
+  // Compute the Gaussian derivatives Lx and Ly
+  image_derivatives_scharr(gaussian,Lx,1,0);
+  image_derivatives_scharr(gaussian,Ly,0,1);
+
+  // Skip the borders for computing the histogram
+  for (int i = 1; i < gaussian.rows-1; i++) {
+    for (int j = 1; j < gaussian.cols-1; j++) {
+      lx = *(Lx.ptr<float>(i)+j);
+      ly = *(Ly.ptr<float>(i)+j);
+      modg = sqrt(lx*lx + ly*ly);
+
+      // Get the maximum
+      if (modg > hmax) {
+        hmax = modg;
+      }
+    }
+  }
+
+  // Skip the borders for computing the histogram
+  for (int i = 1; i < gaussian.rows-1; i++) {
+    for (int j = 1; j < gaussian.cols-1; j++) {
+      lx = *(Lx.ptr<float>(i)+j);
+      ly = *(Ly.ptr<float>(i)+j);
+      modg = sqrt(lx*lx + ly*ly);
+
+      // Find the correspondent bin
+      if (modg != 0.0) {
+        nbin = floor(nbins*(modg/hmax));
+
+        if (nbin == nbins) {
+          nbin--;
+        }
+
+        hist[nbin]++;
+        npoints++;
+      }
+    }
+  }
+
+  // Now find the perc of the histogram percentile
+  nthreshold = (size_t)(npoints*perc);
+
+  for (k = 0; nelements < nthreshold && k < nbins; k++) {
+    nelements = nelements + hist[k];
+  }
+
+  if (nelements < nthreshold) {
+    kperc = 0.03;
+  }
+  else {
+    kperc = hmax*((float)(k)/(float)nbins);
+  }
+
+  delete [] hist;
+  return kperc;
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function computes Scharr image derivatives
+ * @param src Input image
+ * @param dst Output image
+ * @param xorder Derivative order in X-direction (horizontal)
+ * @param yorder Derivative order in Y-direction (vertical)
+ * @param scale Scale factor for the derivative size
+ */
+void compute_scharr_derivatives(const cv::Mat& src, cv::Mat& dst, const size_t& xorder,
+                                const size_t& yorder, const size_t& scale) {
+
+  Mat kx, ky;
+  compute_derivative_kernels(kx, ky, xorder,yorder,scale);
+  sepFilter2D(src,dst,CV_32F,kx,ky);
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function performs a scalar non-linear diffusion step
+ * @param Ld2 Output image in the evolution
+ * @param c Conductivity image
+ * @param Lstep Previous image in the evolution
+ * @param stepsize The step size in time units
+ * @note Forward Euler Scheme 3x3 stencil
+ * The function c is a scalar value that depends on the gradient norm
+ * dL_by_ds = d(c dL_by_dx)_by_dx + d(c dL_by_dy)_by_dy
+ */
+void nld_step_scalar(cv::Mat& Ld, const cv::Mat& c, cv::Mat& Lstep, const float& stepsize) {
+
+#ifdef _OPENMP
+#pragma omp parallel for schedule(dynamic)
+#endif
+  for (int i = 1; i < Lstep.rows-1; i++) {
+    for (int j = 1; j < Lstep.cols-1; j++) {
+      float xpos = ((*(c.ptr<float>(i)+j))+(*(c.ptr<float>(i)+j+1)))*((*(Ld.ptr<float>(i)+j+1))-(*(Ld.ptr<float>(i)+j)));
+      float xneg = ((*(c.ptr<float>(i)+j-1))+(*(c.ptr<float>(i)+j)))*((*(Ld.ptr<float>(i)+j))-(*(Ld.ptr<float>(i)+j-1)));
+
+      float ypos = ((*(c.ptr<float>(i)+j))+(*(c.ptr<float>(i+1)+j)))*((*(Ld.ptr<float>(i+1)+j))-(*(Ld.ptr<float>(i)+j)));
+      float yneg = ((*(c.ptr<float>(i-1)+j))+(*(c.ptr<float>(i)+j)))*((*(Ld.ptr<float>(i)+j))-(*(Ld.ptr<float>(i-1)+j)));
+
+      *(Lstep.ptr<float>(i)+j) = 0.5*stepsize*(xpos-xneg + ypos-yneg);
+    }
+  }
+
+  for (int j = 1; j < Lstep.cols-1; j++) {
+    float xpos = ((*(c.ptr<float>(0)+j))+(*(c.ptr<float>(0)+j+1)))*((*(Ld.ptr<float>(0)+j+1))-(*(Ld.ptr<float>(0)+j)));
+    float xneg = ((*(c.ptr<float>(0)+j-1))+(*(c.ptr<float>(0)+j)))*((*(Ld.ptr<float>(0)+j))-(*(Ld.ptr<float>(0)+j-1)));
+
+    float ypos = ((*(c.ptr<float>(0)+j))+(*(c.ptr<float>(1)+j)))*((*(Ld.ptr<float>(1)+j))-(*(Ld.ptr<float>(0)+j)));
+    float yneg = ((*(c.ptr<float>(0)+j))+(*(c.ptr<float>(0)+j)))*((*(Ld.ptr<float>(0)+j))-(*(Ld.ptr<float>(0)+j)));
+
+    *(Lstep.ptr<float>(0)+j) = 0.5*stepsize*(xpos-xneg + ypos-yneg);
+  }
+
+  for (int j = 1; j < Lstep.cols-1; j++) {
+    float xpos = ((*(c.ptr<float>(Lstep.rows-1)+j))+(*(c.ptr<float>(Lstep.rows-1)+j+1)))*((*(Ld.ptr<float>(Lstep.rows-1)+j+1))-(*(Ld.ptr<float>(Lstep.rows-1)+j)));
+    float xneg = ((*(c.ptr<float>(Lstep.rows-1)+j-1))+(*(c.ptr<float>(Lstep.rows-1)+j)))*((*(Ld.ptr<float>(Lstep.rows-1)+j))-(*(Ld.ptr<float>(Lstep.rows-1)+j-1)));
+
+    float ypos = ((*(c.ptr<float>(Lstep.rows-1)+j))+(*(c.ptr<float>(Lstep.rows-1)+j)))*((*(Ld.ptr<float>(Lstep.rows-1)+j))-(*(Ld.ptr<float>(Lstep.rows-1)+j)));
+    float yneg = ((*(c.ptr<float>(Lstep.rows-2)+j))+(*(c.ptr<float>(Lstep.rows-1)+j)))*((*(Ld.ptr<float>(Lstep.rows-1)+j))-(*(Ld.ptr<float>(Lstep.rows-2)+j)));
+
+    *(Lstep.ptr<float>(Lstep.rows-1)+j) = 0.5*stepsize*(xpos-xneg + ypos-yneg);
+  }
+
+  for (int i = 1; i < Lstep.rows-1; i++) {
+    float xpos = ((*(c.ptr<float>(i)))+(*(c.ptr<float>(i)+1)))*((*(Ld.ptr<float>(i)+1))-(*(Ld.ptr<float>(i))));
+    float xneg = ((*(c.ptr<float>(i)))+(*(c.ptr<float>(i))))*((*(Ld.ptr<float>(i)))-(*(Ld.ptr<float>(i))));
+
+    float ypos = ((*(c.ptr<float>(i)))+(*(c.ptr<float>(i+1))))*((*(Ld.ptr<float>(i+1)))-(*(Ld.ptr<float>(i))));
+    float yneg = ((*(c.ptr<float>(i-1)))+(*(c.ptr<float>(i))))*((*(Ld.ptr<float>(i)))-(*(Ld.ptr<float>(i-1))));
+
+    *(Lstep.ptr<float>(i)) = 0.5*stepsize*(xpos-xneg + ypos-yneg);
+  }
+
+  for (int i = 1; i < Lstep.rows-1; i++) {
+    float xpos = ((*(c.ptr<float>(i)+Lstep.cols-1))+(*(c.ptr<float>(i)+Lstep.cols-1)))*((*(Ld.ptr<float>(i)+Lstep.cols-1))-(*(Ld.ptr<float>(i)+Lstep.cols-1)));
+    float xneg = ((*(c.ptr<float>(i)+Lstep.cols-2))+(*(c.ptr<float>(i)+Lstep.cols-1)))*((*(Ld.ptr<float>(i)+Lstep.cols-1))-(*(Ld.ptr<float>(i)+Lstep.cols-2)));
+
+    float ypos = ((*(c.ptr<float>(i)+Lstep.cols-1))+(*(c.ptr<float>(i+1)+Lstep.cols-1)))*((*(Ld.ptr<float>(i+1)+Lstep.cols-1))-(*(Ld.ptr<float>(i)+Lstep.cols-1)));
+    float yneg = ((*(c.ptr<float>(i-1)+Lstep.cols-1))+(*(c.ptr<float>(i)+Lstep.cols-1)))*((*(Ld.ptr<float>(i)+Lstep.cols-1))-(*(Ld.ptr<float>(i-1)+Lstep.cols-1)));
+
+    *(Lstep.ptr<float>(i)+Lstep.cols-1) = 0.5*stepsize*(xpos-xneg + ypos-yneg);
+  }
+
+  Ld = Ld + Lstep;
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function downsamples the input image with the kernel [1/4,1/2,1/4]
+ * @param img Input image to be downsampled
+ * @param dst Output image with half of the resolution of the input image
+ */
+void downsample_image(const cv::Mat& src, cv::Mat& dst) {
+
+  int i1 = 0, j1 = 0, i2 = 0, j2 = 0;
+
+  for (i1 = 1; i1 < src.rows; i1+=2) {
+    j2 = 0;
+    for (j1 = 1; j1 < src.cols; j1+=2) {
+      *(dst.ptr<float>(i2)+j2) = 0.5*(*(src.ptr<float>(i1)+j1))+0.25*(*(src.ptr<float>(i1)+j1-1) + *(src.ptr<float>(i1)+j1+1));
+      j2++;
+    }
+
+    i2++;
+  }
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function downsamples the input image using OpenCV resize
+ * @param img Input image to be downsampled
+ * @param dst Output image with half of the resolution of the input image
+ */
+void halfsample_image(const cv::Mat& src, cv::Mat& dst) {
+
+  // Make sure the destination image is of the right size
+  CV_Assert(src.cols/2==dst.cols);
+  CV_Assert(src.rows / 2 == dst.rows);
+  resize(src,dst,dst.size(),0,0,cv::INTER_AREA);
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief Compute Scharr derivative kernels for sizes different than 3
+ * @param kx_ The derivative kernel in x-direction
+ * @param ky_ The derivative kernel in y-direction
+ * @param dx The derivative order in x-direction
+ * @param dy The derivative order in y-direction
+ * @param scale The kernel size
+ */
+void compute_derivative_kernels(cv::OutputArray kx_, cv::OutputArray ky_,
+                                const size_t& dx, const size_t& dy, const size_t& scale) {
+
+  const int ksize = 3 + 2*(scale-1);
+
+  // The usual Scharr kernel
+  if (scale == 1) {
+    getDerivKernels(kx_,ky_,dx,dy,0,true,CV_32F);
+    return;
+  }
+
+  kx_.create(ksize,1,CV_32F,-1,true);
+  ky_.create(ksize,1,CV_32F,-1,true);
+  Mat kx = kx_.getMat();
+  Mat ky = ky_.getMat();
+
+  float w = 10.0/3.0;
+  float norm = 1.0/(2.0*scale*(w+2.0));
+
+  for (int k = 0; k < 2; k++) {
+    Mat* kernel = k == 0 ? &kx : &ky;
+    int order = k == 0 ? dx : dy;
+    float kerI[1000];
+
+    for (int t = 0; t<ksize; t++) {
+      kerI[t] = 0;
+    }
+
+    if (order == 0) {
+      kerI[0] = norm;
+      kerI[ksize/2] = w*norm;
+      kerI[ksize-1] = norm;
+    }
+    else if (order == 1) {
+      kerI[0] = -1;
+      kerI[ksize/2] = 0;
+      kerI[ksize-1] = 1;
+    }
+
+    Mat temp(kernel->rows, kernel->cols, CV_32F, &kerI[0]);
+    temp.copyTo(*kernel);
+  }
+}

modules/features2d/src/akaze/nldiffusion_functions.h

@@ -0,0 +1,41 @@
+#ifndef _NLDIFFUSION_FUNCTIONS_H_
+#define _NLDIFFUSION_FUNCTIONS_H_
+
+//******************************************************************************
+//******************************************************************************
+
+// Includes
+#include "precomp.hpp"
+
+// OpenMP Includes
+#ifdef _OPENMP
+# include <omp.h>
+#endif
+
+//*************************************************************************************
+//*************************************************************************************
+
+// Declaration of functions
+void gaussian_2D_convolution(const cv::Mat& src, cv::Mat& dst, const size_t& ksize_x,
+                             const size_t& ksize_y, const float& sigma);
+void image_derivatives_scharr(const cv::Mat& src, cv::Mat& dst,
+                              const size_t& xorder, const size_t& yorder);
+void pm_g1(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k);
+void pm_g2(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k);
+void weickert_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k);
+void charbonnier_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k);
+float compute_k_percentile(const cv::Mat& img, const float& perc, const float& gscale,
+                           const size_t& nbins, const size_t& ksize_x, const size_t& ksize_y);
+void compute_scharr_derivatives(const cv::Mat& src, cv::Mat& dst, const size_t& xorder,
+                                const size_t& yorder, const size_t& scale);
+void nld_step_scalar(cv::Mat& Ld, const cv::Mat& c, cv::Mat& Lstep, const float& stepsize);
+void downsample_image(const cv::Mat& src, cv::Mat& dst);
+void halfsample_image(const cv::Mat& src, cv::Mat& dst);
+void compute_derivative_kernels(cv::OutputArray kx_, cv::OutputArray ky_,
+                                const size_t& dx, const size_t& dy, const size_t& scale);
+
+//*************************************************************************************
+//*************************************************************************************
+
+
+#endif

modules/features2d/src/akaze/utils.cpp

@@ -0,0 +1,196 @@
+//=============================================================================
+//
+// utils.cpp
+// Authors: Pablo F. Alcantarilla (1), Jesus Nuevo (2)
+// Institutions: Georgia Institute of Technology (1)
+//               TrueVision Solutions (2)
+//
+// Date: 15/09/2013
+// Email: pablofdezalc@gmail.com
+//
+// AKAZE Features Copyright 2013, Pablo F. Alcantarilla, Jesus Nuevo
+// All Rights Reserved
+// See LICENSE for the license information
+//=============================================================================
+
+/**
+ * @file utils.cpp
+ * @brief Some utilities functions
+ * @date Sep 15, 2013
+ * @author Pablo F. Alcantarilla, Jesus Nuevo
+ */
+
+#include "precomp.hpp"
+#include "utils.h"
+
+// Namespaces
+using namespace std;
+using namespace cv;
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function computes the minimum value of a float image
+ * @param src Input image
+ * @param value Minimum value
+ */
+void compute_min_32F(const cv::Mat &src, float &value) {
+
+  float aux = 1000.0;
+
+  for (int i = 0; i < src.rows; i++) {
+    for (int j = 0; j < src.cols; j++) {
+      if (src.at<float>(i,j) < aux) {
+        aux = src.at<float>(i,j);
+      }
+    }
+  }
+
+  value = aux;
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function computes the maximum value of a float image
+ * @param src Input image
+ * @param value Maximum value
+ */
+void compute_max_32F(const cv::Mat &src, float &value) {
+
+  float aux = 0.0;
+
+  for (int i = 0; i < src.rows; i++) {
+    for (int j = 0; j < src.cols; j++) {
+      if (src.at<float>(i,j) > aux) {
+        aux = src.at<float>(i,j);
+      }
+    }
+  }
+
+  value = aux;
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function converts the scale of the input image prior to visualization
+ * @param src Input/Output image
+ * @param value Maximum value
+ */
+void convert_scale(cv::Mat &src) {
+
+  float min_val = 0, max_val = 0;
+
+  compute_min_32F(src,min_val);
+
+  src = src - min_val;
+
+  compute_max_32F(src,max_val);
+  src = src / max_val;
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function copies the input image and converts the scale of the copied
+ * image prior visualization
+ * @param src Input image
+ * @param dst Output image
+ */
+void copy_and_convert_scale(const cv::Mat &src, cv::Mat dst) {
+
+  float min_val = 0, max_val = 0;
+
+  src.copyTo(dst);
+  compute_min_32F(dst,min_val);
+
+  dst = dst - min_val;
+
+  compute_max_32F(dst,max_val);
+  dst = dst / max_val;
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+const size_t length = string("--descriptor_channels").size() + 2;
+static inline std::ostream& cout_help()
+{ cout << setw(length); return cout; }
+
+static inline std::string toUpper(std::string s)
+{
+  std::transform(s.begin(), s.end(), s.begin(), ::toupper);
+  return s;
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function shows the possible command line configuration options
+ */
+void show_input_options_help(int example) {
+
+  fflush(stdout);
+  cout << "A-KAZE Features" << endl;
+  cout << "Usage: ";
+  if (example == 0) {
+    cout << "./akaze_features -i img.jpg [options]" << endl;
+  }
+  else if (example == 1) {
+    cout << "./akaze_match img1.jpg img2.pgm homography.txt [options]" << endl;
+  }
+  else if (example == 2) {
+    cout << "./akaze_compare img1.jpg img2.pgm homography.txt [options]" << endl;
+  }
+  
+  cout << endl;
+  cout_help() << "Options below are not mandatory. Unless specified, default arguments are used." << endl << endl;  
+  // Justify on the left
+  cout << left;
+  // Generalities
+  cout_help() << "--help" << "Show the command line options" << endl;
+  cout_help() << "--verbose " << "Verbosity is required" << endl;
+  cout_help() << endl;
+  // Scale-space parameters
+  cout_help() << "--soffset" << "Base scale offset (sigma units)" << endl;
+  cout_help() << "--omax" << "Maximum octave of image evolution" << endl;
+  cout_help() << "--nsublevels" << "Number of sublevels per octave" << endl;
+  cout_help() << "--diffusivity" << "Diffusivity function. Possible values:" << endl;
+  cout_help() << " " << "0 -> Perona-Malik, g1 = exp(-|dL|^2/k^2)" << endl;
+  cout_help() << " " << "1 -> Perona-Malik, g2 = 1 / (1 + dL^2 / k^2)" << endl;
+  cout_help() << " " << "2 -> Weickert diffusivity" << endl;
+  cout_help() << " " << "3 -> Charbonnier diffusivity" << endl;
+  cout_help() << endl;
+  // Feature detection parameters.
+  cout_help() << "--dthreshold" << "Feature detector threshold response for keypoints" << endl;
+  cout_help() << " " << "(0.001 can be a good value)" << endl;
+  cout_help() << endl;
+  // Descriptor parameters.
+  cout_help() << "--descriptor" << "Descriptor Type. Possible values:" << endl;
+  cout_help() << " " << "0 -> SURF_UPRIGHT" << endl;
+  cout_help() << " " << "1 -> SURF" << endl;
+  cout_help() << " " << "2 -> M-SURF_UPRIGHT," << endl;
+  cout_help() << " " << "3 -> M-SURF" << endl;
+  cout_help() << " " << "4 -> M-LDB_UPRIGHT" << endl;
+  cout_help() << " " << "5 -> M-LDB" << endl;
+  
+  cout_help() << "--descriptor_channels " << "Descriptor Channels for M-LDB. Valid values: " << endl;
+  cout_help() << " " << "1 -> intensity" << endl;
+  cout_help() << " " << "2 -> intensity + gradient magnitude" << endl;
+  cout_help() << " " << "3 -> intensity + X and Y gradients" <<endl;
+  
+  cout_help() << "--descriptor_size" << "Descriptor size for M-LDB in bits." << endl;
+  cout_help() << " " << "0: means the full length descriptor (486)!!" << endl;
+  cout_help() << endl;
+  // Save results?
+  cout_help() << "--show_results" << "Possible values below:" << endl;
+  cout_help() << " " << "1 -> show detection results." << endl;
+  cout_help() << " " << "0 -> don't show detection results" << endl;
+  cout_help() << endl;
+}

modules/features2d/src/akaze/utils.h

@@ -0,0 +1,54 @@
+
+#ifndef _UTILS_H_
+#define _UTILS_H_
+
+//******************************************************************************
+//******************************************************************************
+
+// OpenCV Includes
+#include "precomp.hpp"
+
+// System Includes
+#include <stdlib.h>
+#include <stdio.h>
+#include <cstdlib>
+#include <vector>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+
+//******************************************************************************
+//******************************************************************************
+
+// Stringify common types such as int, double and others.
+template <typename T>
+inline std::string to_string(const T& x) {
+  std::stringstream oss;
+  oss << x;
+  return oss.str();
+}
+
+//******************************************************************************
+//******************************************************************************
+
+// Stringify and format integral types as follows:
+// to_formatted_string(  1, 2) produces string:  '01'
+// to_formatted_string(  5, 2) produces string:  '05'
+// to_formatted_string( 19, 2) produces string:  '19'
+// to_formatted_string( 19, 3) produces string: '019'
+template <typename Integer>
+inline std::string to_formatted_string(Integer x, int num_digits) {
+  std::stringstream oss;
+  oss << std::setfill('0') << std::setw(num_digits) << x;
+  return oss.str();
+}
+
+//******************************************************************************
+//******************************************************************************
+
+void compute_min_32F(const cv::Mat& src, float& value);
+void compute_max_32F(const cv::Mat& src, float& value);
+void convert_scale(cv::Mat& src);
+void copy_and_convert_scale(const cv::Mat& src, cv::Mat& dst);
+
+#endif

modules/features2d/src/kaze/KAZE.h

@@ -0,0 +1,294 @@
+
+/**
+ * @file KAZE.h
+ * @brief Main program for detecting and computing descriptors in a nonlinear
+ * scale space
+ * @date Jan 21, 2012
+ * @author Pablo F. Alcantarilla
+ */
+
+#ifndef KAZE_H_
+#define KAZE_H_
+
+//*************************************************************************************
+//*************************************************************************************
+
+// Includes
+#include "config.h"
+#include "nldiffusion_functions.h"
+#include "fed.h"
+#include "utils.h"
+
+//*************************************************************************************
+//*************************************************************************************
+
+// KAZE Class Declaration
+class KAZE {
+
+private:
+
+  // Parameters of the Nonlinear diffusion class
+  float soffset_;	// Base scale offset
+  float sderivatives_; // Standard deviation of the Gaussian for the nonlinear diff. derivatives
+  int omax_;		// Maximum octave level
+  int nsublevels_;	// Number of sublevels per octave level
+  int img_width_;	// Width of the original image
+  int img_height_; // Height of the original image
+  bool save_scale_space_; // For saving scale space images
+  bool verbosity_;	// Verbosity level
+  std::vector<TEvolution> evolution_;	// Vector of nonlinear diffusion evolution
+  float kcontrast_; // The contrast parameter for the scalar nonlinear diffusion
+  float dthreshold_;	// Feature detector threshold response
+  int diffusivity_; 	// Diffusivity type, 0->PM G1, 1->PM G2, 2-> Weickert
+  int descriptor_mode_; // Descriptor mode
+  bool use_fed_;        // Set to true in case we want to use FED for the nonlinear diffusion filtering. Set false for using AOS
+  bool use_upright_;	// Set to true in case we want to use the upright version of the descriptors
+  bool use_extended_;	// Set to true in case we want to use the extended version of the descriptors
+
+  // Vector of keypoint vectors for finding extrema in multiple threads
+  std::vector<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
+
+  // Computation times variables in ms
+  double tkcontrast_;       // Kcontrast factor computation
+  double tnlscale_;         // Nonlinear Scale space generation
+  double tdetector_;        // Feature detector
+  double tmderivatives_;    // Multiscale derivatives computation
+  double tdresponse_;       // Detector response computation
+  double tdescriptor_;      // Feature descriptor
+  double tsubpixel_;        // Subpixel refinement
+
+  // Some auxiliary variables used in the AOS step
+  cv::Mat Ltx_, Lty_, px_, py_, ax_, ay_, bx_, by_, qr_, qc_;
+
+public:
+
+  // Constructor
+  KAZE(KAZEOptions& options);
+
+  // Destructor
+  ~KAZE(void);
+
+  // Public methods for KAZE interface
+  void Allocate_Memory_Evolution(void);
+  int Create_Nonlinear_Scale_Space(const cv::Mat& img);
+  void Feature_Detection(std::vector<cv::KeyPoint>& kpts);
+  void Feature_Description(std::vector<cv::KeyPoint>& kpts, cv::Mat& desc);
+
+  // Methods for saving the scale space set of images and detector responses
+  void Save_Nonlinear_Scale_Space(void);
+  void Save_Detector_Responses(void);
+  void Save_Flow_Responses(void);
+
+private:
+
+  // Feature Detection Methods
+  void Compute_KContrast(const cv::Mat& img, const float& kper);
+  void Compute_Multiscale_Derivatives(void);
+  void Compute_Detector_Response(void);
+  void Determinant_Hessian_Parallel(std::vector<cv::KeyPoint>& kpts);
+  void Find_Extremum_Threading(const int& level);
+  void Do_Subpixel_Refinement(std::vector<cv::KeyPoint>& kpts);
+  void Feature_Suppression_Distance(std::vector<cv::KeyPoint>& kpts, const float& mdist);
+
+  // AOS Methods
+  void AOS_Step_Scalar(cv::Mat &Ld, const cv::Mat &Ldprev, const cv::Mat &c, const float& stepsize);
+  void AOS_Rows(const cv::Mat &Ldprev, const cv::Mat &c, const float& stepsize);
+  void AOS_Columns(const cv::Mat &Ldprev, const cv::Mat &c, const float& stepsize);
+  void Thomas(const cv::Mat &a, const cv::Mat &b, const cv::Mat &Ld, cv::Mat &x);
+
+  // Feature Description methods
+  void Compute_Main_Orientation_SURF(cv::KeyPoint& kpt);
+
+  // Descriptor Mode -> 0 SURF 64
+  void Get_SURF_Upright_Descriptor_64(const cv::KeyPoint& kpt, float* desc);
+  void Get_SURF_Descriptor_64(const cv::KeyPoint& kpt, float* desc);
+
+  // Descriptor Mode -> 0 SURF 128
+  void Get_SURF_Upright_Descriptor_128(const cv::KeyPoint& kpt, float* desc);
+  void Get_SURF_Descriptor_128(const cv::KeyPoint& kpt, float* desc);
+
+  // Descriptor Mode -> 1 M-SURF 64
+  void Get_MSURF_Upright_Descriptor_64(const cv::KeyPoint& kpt, float* desc);
+  void Get_MSURF_Descriptor_64(const cv::KeyPoint& kpt, float* desc);
+
+  // Descriptor Mode -> 1 M-SURF 128
+  void Get_MSURF_Upright_Descriptor_128(const cv::KeyPoint& kpt, float* desc);
+  void Get_MSURF_Descriptor_128(const cv::KeyPoint& kpt, float *desc);
+
+  // Descriptor Mode -> 2 G-SURF 64
+  void Get_GSURF_Upright_Descriptor_64(const cv::KeyPoint& kpt, float* desc);
+  void Get_GSURF_Descriptor_64(const cv::KeyPoint& kpt, float *desc);
+
+  // Descriptor Mode -> 2 G-SURF 128
+  void Get_GSURF_Upright_Descriptor_128(const cv::KeyPoint& kpt, float* desc);
+  void Get_GSURF_Descriptor_128(const cv::KeyPoint& kpt, float* desc);
+
+public:
+
+  // Setters
+  void Set_Scale_Offset(float soffset) {
+    soffset_ = soffset;
+  }
+
+  void Set_SDerivatives(float sderivatives) {
+    sderivatives_ = sderivatives;
+  }
+
+  void Set_Octave_Max(int omax) {
+    omax_ = omax;
+  }
+
+  void Set_NSublevels(int nsublevels) {
+    nsublevels_ = nsublevels;
+  }
+
+  void Set_Save_Scale_Space_Flag(bool save_scale_space) {
+    save_scale_space_ = save_scale_space;
+  }
+
+  void Set_Image_Width(int img_width) {
+    img_width_ = img_width;
+  }
+
+  void Set_Image_Height(int img_height) {
+    img_height_ = img_height;
+  }
+
+  void Set_Verbosity_Level(bool verbosity) {
+    verbosity_ = verbosity;
+  }
+
+  void Set_KContrast(float kcontrast) {
+    kcontrast_ = kcontrast;
+  }
+
+  void Set_Detector_Threshold(float dthreshold) {
+    dthreshold_ = dthreshold;
+  }
+
+  void Set_Diffusivity_Type(int diffusivity) {
+    diffusivity_ = diffusivity;
+  }
+
+  void Set_Descriptor_Mode(int descriptor_mode) {
+    descriptor_mode_ = descriptor_mode;
+  }
+
+  void Set_Use_FED(bool use_fed) {
+    use_fed_ = use_fed;
+  }
+
+  void Set_Upright(bool use_upright) {
+    use_upright_ = use_upright;
+  }
+
+  void Set_Extended(bool use_extended) {
+    use_extended_ = use_extended;
+  }
+
+  // Getters
+  float Get_Scale_Offset(void) {
+    return soffset_;
+  }
+
+  float Get_SDerivatives(void) {
+    return sderivatives_;
+  }
+
+  int Get_Octave_Max(void) {
+    return omax_;
+  }
+
+  int Get_NSublevels(void) {
+    return nsublevels_;
+  }
+
+  bool Get_Save_Scale_Space_Flag(void) {
+    return save_scale_space_;
+  }
+
+  int Get_Image_Width(void) {
+    return img_width_;
+  }
+
+  int Get_Image_Height(void) {
+    return img_height_;
+  }
+
+  bool Get_Verbosity_Level(void) {
+    return verbosity_;
+  }
+
+  float Get_KContrast(void) {
+    return kcontrast_;
+  }
+
+  float Get_Detector_Threshold(void) {
+    return dthreshold_;
+  }
+
+  int Get_Diffusivity_Type(void) {
+    return diffusivity_;
+  }
+
+  int Get_Descriptor_Mode(void) {
+    return descriptor_mode_;
+  }
+
+  bool Get_Upright(void) {
+    return use_upright_;
+  }
+
+  bool Get_Extended(void) {
+    return use_extended_;
+  }
+
+  float Get_Time_KContrast(void) {
+    return tkcontrast_;
+  }
+
+  float Get_Time_NLScale(void) {
+    return tnlscale_;
+  }
+
+  float Get_Time_Detector(void) {
+    return tdetector_;
+  }
+
+  float Get_Time_Multiscale_Derivatives(void) {
+    return tmderivatives_;
+  }
+
+  float Get_Time_Detector_Response(void) {
+    return tdresponse_;
+  }
+
+  float Get_Time_Descriptor(void) {
+    return tdescriptor_;
+  }
+
+  float Get_Time_Subpixel(void) {
+    return tsubpixel_;
+  }
+};
+
+//*************************************************************************************
+//*************************************************************************************
+
+// Inline functions
+float getAngle(const float& x, const float& y);
+float gaussian(const float& x, const float& y, const float& sig);
+void checkDescriptorLimits(int &x, int &y, const int& width, const int& height);
+void clippingDescriptor(float *desc, const int& dsize, const int& niter, const float& ratio);
+int fRound(const float& flt);
+
+//*************************************************************************************
+//*************************************************************************************
+
+#endif // KAZE_H_

modules/features2d/src/kaze/config.h

@@ -0,0 +1,129 @@
+
+/**
+ * @file config.h
+ * @brief Configuration file
+ * @date Dec 27, 2011
+ * @author Pablo F. Alcantarilla
+ */
+
+#ifndef _CONFIG_H_
+#define _CONFIG_H_
+
+//******************************************************************************
+//******************************************************************************
+
+// System Includes
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <cstdlib>
+#include <string>
+#include <vector>
+#include <math.h>
+
+// OpenCV Includes
+#include "precomp.hpp"
+
+// OpenMP Includes
+#ifdef _OPENMP
+#include <omp.h>
+#else
+#define omp_get_thread_num() 0
+#endif
+
+//*************************************************************************************
+//*************************************************************************************
+
+// Some defines
+#define NMAX_CHAR 400
+
+// Some default options
+const float DEFAULT_SCALE_OFFSET = 1.60; // Base scale offset (sigma units)
+const float DEFAULT_OCTAVE_MAX = 4.0; // Maximum octave evolution of the image 2^sigma (coarsest scale sigma units)
+const int DEFAULT_NSUBLEVELS = 4; // Default number of sublevels per scale level
+const float DEFAULT_DETECTOR_THRESHOLD = 0.001; // Detector response threshold to accept point
+const float DEFAULT_MIN_DETECTOR_THRESHOLD = 0.00001;     // Minimum Detector response threshold to accept point
+const int DEFAULT_DESCRIPTOR_MODE = 1; // Descriptor Mode 0->SURF, 1->M-SURF
+const bool DEFAULT_USE_FED = true;  // 0->AOS, 1->FED
+const bool DEFAULT_UPRIGHT = false;  // Upright descriptors, not invariant to rotation
+const bool DEFAULT_EXTENDED = false; // Extended descriptor, dimension 128
+const bool DEFAULT_SAVE_SCALE_SPACE = false; // For saving the scale space images
+const bool DEFAULT_VERBOSITY = false; // Verbosity level (0->no verbosity)
+const bool DEFAULT_SHOW_RESULTS = true; // For showing the output image with the detected features plus some ratios
+const bool DEFAULT_SAVE_KEYPOINTS = false; // For saving the list of keypoints
+
+// Some important configuration variables
+const float DEFAULT_SIGMA_SMOOTHING_DERIVATIVES = 1.0;
+const float DEFAULT_KCONTRAST = .01;
+const float KCONTRAST_PERCENTILE = 0.7;
+const int KCONTRAST_NBINS = 300;
+const bool COMPUTE_KCONTRAST = true;
+const int DEFAULT_DIFFUSIVITY_TYPE = 1;  // 0 -> PM G1, 1 -> PM G2, 2 -> Weickert
+const bool USE_CLIPPING_NORMALIZATION = false;
+const float CLIPPING_NORMALIZATION_RATIO = 1.6;
+const int CLIPPING_NORMALIZATION_NITER = 5;
+
+//*************************************************************************************
+//*************************************************************************************
+
+struct KAZEOptions {
+
+  KAZEOptions() {
+    // Load the default options
+    soffset = DEFAULT_SCALE_OFFSET;
+    omax = DEFAULT_OCTAVE_MAX;
+    nsublevels = DEFAULT_NSUBLEVELS;
+    dthreshold = DEFAULT_DETECTOR_THRESHOLD;
+    use_fed = DEFAULT_USE_FED;
+    upright = DEFAULT_UPRIGHT;
+    extended = DEFAULT_EXTENDED;
+    descriptor = DEFAULT_DESCRIPTOR_MODE;
+    diffusivity = DEFAULT_DIFFUSIVITY_TYPE;
+    sderivatives = DEFAULT_SIGMA_SMOOTHING_DERIVATIVES;
+    save_scale_space = DEFAULT_SAVE_SCALE_SPACE;
+    save_keypoints = DEFAULT_SAVE_KEYPOINTS;
+    verbosity = DEFAULT_VERBOSITY;
+    show_results = DEFAULT_SHOW_RESULTS;
+  }
+
+  float soffset;
+  int omax;
+  int nsublevels;
+  int img_width;
+  int img_height;
+  int diffusivity;
+  float sderivatives;
+  float dthreshold;
+  bool use_fed;
+  bool upright;
+  bool extended;
+  int descriptor;
+  bool save_scale_space;
+  bool save_keypoints;
+  bool verbosity;
+  bool show_results;
+};
+
+struct TEvolution {
+  cv::Mat Lx, Ly;	// First order spatial derivatives
+  cv::Mat Lxx, Lxy, Lyy;	// Second order spatial derivatives
+  cv::Mat Lflow;	// Diffusivity image
+  cv::Mat Lt;	// Evolution image
+  cv::Mat Lsmooth; // Smoothed image
+  cv::Mat Lstep; // Evolution step update
+  cv::Mat Ldet; // Detector response
+  float etime;	// Evolution time
+  float esigma;	// Evolution sigma. For linear diffusion t = sigma^2 / 2
+  float octave;	// Image octave
+  float sublevel;	// Image sublevel in each octave
+  int sigma_size;	// Integer esigma. For computing the feature detector responses
+};
+
+//*************************************************************************************
+//*************************************************************************************
+
+#endif
+
+
+
+

modules/features2d/src/kaze/fed.cpp

@@ -0,0 +1,192 @@
+//=============================================================================
+//
+// fed.cpp
+// Authors: Pablo F. Alcantarilla (1), Jesus Nuevo (2)
+// Institutions: Georgia Institute of Technology (1)
+//               TrueVision Solutions (2)
+// Date: 15/09/2013
+// Email: pablofdezalc@gmail.com
+//
+// AKAZE Features Copyright 2013, Pablo F. Alcantarilla, Jesus Nuevo
+// All Rights Reserved
+// See LICENSE for the license information
+//=============================================================================
+
+/**
+ * @file fed.cpp
+ * @brief Functions for performing Fast Explicit Diffusion and building the
+ * nonlinear scale space
+ * @date Sep 15, 2013
+ * @author Pablo F. Alcantarilla, Jesus Nuevo
+ * @note This code is derived from FED/FJ library from Grewenig et al.,
+ * The FED/FJ library allows solving more advanced problems
+ * Please look at the following papers for more information about FED:
+ * [1] S. Grewenig, J. Weickert, C. Schroers, A. Bruhn. Cyclic Schemes for
+ * PDE-Based Image Analysis. Technical Report No. 327, Department of Mathematics,
+ * Saarland University, Saarbrücken, Germany, March 2013
+ * [2] S. Grewenig, J. Weickert, A. Bruhn. From box filtering to fast explicit diffusion.
+ * DAGM, 2010
+ *
+*/
+
+#include "fed.h"
+
+using namespace std;
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function allocates an array of the least number of time steps such
+ * that a certain stopping time for the whole process can be obtained and fills
+ * it with the respective FED time step sizes for one cycle
+ * The function returns the number of time steps per cycle or 0 on failure
+ * @param T Desired process stopping time
+ * @param M Desired number of cycles
+ * @param tau_max Stability limit for the explicit scheme
+ * @param reordering Reordering flag
+ * @param tau The vector with the dynamic step sizes
+ */
+int fed_tau_by_process_time(const float& T, const int& M, const float& tau_max,
+                            const bool& reordering, std::vector<float>& tau) {
+  // All cycles have the same fraction of the stopping time
+  return fed_tau_by_cycle_time(T/(float)M,tau_max,reordering,tau);
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function allocates an array of the least number of time steps such
+ * that a certain stopping time for the whole process can be obtained and fills it
+ * it with the respective FED time step sizes for one cycle
+ * The function returns the number of time steps per cycle or 0 on failure
+ * @param t Desired cycle stopping time
+ * @param tau_max Stability limit for the explicit scheme
+ * @param reordering Reordering flag
+ * @param tau The vector with the dynamic step sizes
+ */
+int fed_tau_by_cycle_time(const float& t, const float& tau_max,
+                          const bool& reordering, std::vector<float> &tau) {
+  int n = 0;          // Number of time steps
+  float scale = 0.0;  // Ratio of t we search to maximal t
+
+  // Compute necessary number of time steps
+  n = (int)(ceilf(sqrtf(3.0*t/tau_max+0.25f)-0.5f-1.0e-8f)+ 0.5f);
+  scale = 3.0*t/(tau_max*(float)(n*(n+1)));
+
+  // Call internal FED time step creation routine
+  return fed_tau_internal(n,scale,tau_max,reordering,tau);
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function allocates an array of time steps and fills it with FED
+ * time step sizes
+ * The function returns the number of time steps per cycle or 0 on failure
+ * @param n Number of internal steps
+ * @param scale Ratio of t we search to maximal t
+ * @param tau_max Stability limit for the explicit scheme
+ * @param reordering Reordering flag
+ * @param tau The vector with the dynamic step sizes
+ */
+int fed_tau_internal(const int& n, const float& scale, const float& tau_max,
+                     const bool& reordering, std::vector<float> &tau) {
+  float c = 0.0, d = 0.0;     // Time savers
+  vector<float> tauh;    // Helper vector for unsorted taus
+
+  if (n <= 0) {
+    return 0;
+  }
+
+  // Allocate memory for the time step size
+  tau = vector<float>(n);
+
+  if (reordering) {
+    tauh = vector<float>(n);
+  }
+
+  // Compute time saver
+  c = 1.0f / (4.0f * (float)n + 2.0f);
+  d = scale * tau_max / 2.0f;
+
+  // Set up originally ordered tau vector
+  for (int k = 0; k < n; ++k) {
+    float h = cosf(CV_PI * (2.0f * (float)k + 1.0f) * c);
+
+    if (reordering) {
+      tauh[k] = d / (h * h);
+    }
+    else {
+      tau[k] = d / (h * h);
+    }
+  }
+
+  // Permute list of time steps according to chosen reordering function
+  int kappa = 0, prime = 0;
+
+  if (reordering == true) {
+    // Choose kappa cycle with k = n/2
+    // This is a heuristic. We can use Leja ordering instead!!
+    kappa = n / 2;
+
+    // Get modulus for permutation
+    prime = n + 1;
+
+    while (!fed_is_prime_internal(prime)) {
+      prime++;
+    }
+
+    // Perform permutation
+    for (int k = 0, l = 0; l < n; ++k, ++l) {
+      int index = 0;
+      while ((index = ((k+1)*kappa) % prime - 1) >= n) {
+        k++;
+      }
+
+      tau[l] = tauh[index];
+    }
+  }
+
+  return n;
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function checks if a number is prime or not
+ * @param number Number to check if it is prime or not
+ * @return true if the number is prime
+ */
+bool fed_is_prime_internal(const int& number) {
+  bool is_prime = false;
+
+  if (number <= 1) {
+    return false;
+  }
+  else if (number == 1 || number == 2 || number == 3 || number == 5 || number == 7) {
+    return true;
+  }
+  else if ((number % 2) == 0 || (number % 3) == 0 || (number % 5) == 0 || (number % 7) == 0) {
+    return false;
+  }
+  else {
+    is_prime = true;
+    int upperLimit = sqrt(number+1.0);
+    int divisor = 11;
+
+    while (divisor <= upperLimit ) {
+      if (number % divisor == 0)
+      {
+        is_prime = false;
+      }
+
+      divisor +=2;
+    }
+
+    return is_prime;
+  }
+}

modules/features2d/src/kaze/fed.h

@@ -0,0 +1,30 @@
+#ifndef FED_H
+#define FED_H
+
+//******************************************************************************
+//******************************************************************************
+
+// Includes
+#include <iostream>
+#include <stdlib.h>
+#include <stdio.h>
+#include <cstdlib>
+#include <math.h>
+#include <vector>
+
+//*************************************************************************************
+//*************************************************************************************
+
+// Declaration of functions
+int fed_tau_by_process_time(const float& T, const int& M, const float& tau_max,
+                            const bool& reordering, std::vector<float>& tau);
+int fed_tau_by_cycle_time(const float& t, const float& tau_max,
+                          const bool& reordering, std::vector<float> &tau) ;
+int fed_tau_internal(const int& n, const float& scale, const float& tau_max,
+                     const bool& reordering, std::vector<float> &tau);
+bool fed_is_prime_internal(const int& number);
+
+//*************************************************************************************
+//*************************************************************************************
+
+#endif // FED_H

modules/features2d/src/kaze/nldiffusion_functions.cpp

@@ -0,0 +1,386 @@
+
+//=============================================================================
+//
+// nldiffusion_functions.cpp
+// Author: Pablo F. Alcantarilla
+// Institution: University d'Auvergne
+// Address: Clermont Ferrand, France
+// Date: 27/12/2011
+// Email: pablofdezalc@gmail.com
+//
+// KAZE Features Copyright 2012, Pablo F. Alcantarilla
+// All Rights Reserved
+// See LICENSE for the license information
+//=============================================================================
+
+/**
+ * @file nldiffusion_functions.cpp
+ * @brief Functions for non-linear diffusion applications:
+ * 2D Gaussian Derivatives
+ * Perona and Malik conductivity equations
+ * Perona and Malik evolution
+ * @date Dec 27, 2011
+ * @author Pablo F. Alcantarilla
+ */
+
+#include "nldiffusion_functions.h"
+
+// Namespaces
+using namespace std;
+using namespace cv;
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function smoothes an image with a Gaussian kernel
+ * @param src Input image
+ * @param dst Output image
+ * @param ksize_x Kernel size in X-direction (horizontal)
+ * @param ksize_y Kernel size in Y-direction (vertical)
+ * @param sigma Kernel standard deviation
+ */
+void gaussian_2D_convolution(const cv::Mat& src, cv::Mat& dst,
+                             int ksize_x, int ksize_y, float sigma) {
+
+  size_t ksize_x_ = 0, ksize_y_ = 0;
+
+  // Compute an appropriate kernel size according to the specified sigma
+  if (sigma > ksize_x || sigma > ksize_y || ksize_x == 0 || ksize_y == 0) {
+    ksize_x_ = ceil(2.0*(1.0 + (sigma-0.8)/(0.3)));
+    ksize_y_ = ksize_x_;
+  }
+
+  // The kernel size must be and odd number
+  if ((ksize_x_ % 2) == 0) {
+    ksize_x_ += 1;
+  }
+
+  if ((ksize_y_ % 2) == 0) {
+    ksize_y_ += 1;
+  }
+
+  // Perform the Gaussian Smoothing with border replication
+  GaussianBlur(src,dst,Size(ksize_x_,ksize_y_),sigma,sigma,cv::BORDER_REPLICATE);
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function computes the Perona and Malik conductivity coefficient g1
+ * g1 = exp(-|dL|^2/k^2)
+ * @param Lx First order image derivative in X-direction (horizontal)
+ * @param Ly First order image derivative in Y-direction (vertical)
+ * @param dst Output image
+ * @param k Contrast factor parameter
+ */
+void pm_g1(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k) {
+  cv::exp(-(Lx.mul(Lx) + Ly.mul(Ly))/(k*k),dst);
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function computes the Perona and Malik conductivity coefficient g2
+ * g2 = 1 / (1 + dL^2 / k^2)
+ * @param Lx First order image derivative in X-direction (horizontal)
+ * @param Ly First order image derivative in Y-direction (vertical)
+ * @param dst Output image
+ * @param k Contrast factor parameter
+ */
+void pm_g2(const cv::Mat &Lx, const cv::Mat& Ly, cv::Mat& dst, float k) {
+  dst = 1./(1. + (Lx.mul(Lx) + Ly.mul(Ly))/(k*k));
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function computes Weickert conductivity coefficient g3
+ * @param Lx First order image derivative in X-direction (horizontal)
+ * @param Ly First order image derivative in Y-direction (vertical)
+ * @param dst Output image
+ * @param k Contrast factor parameter
+ * @note For more information check the following paper: J. Weickert
+ * Applications of nonlinear diffusion in image processing and computer vision,
+ * Proceedings of Algorithmy 2000
+ */
+void weickert_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k) {
+  Mat modg;
+  cv::pow((Lx.mul(Lx) + Ly.mul(Ly))/(k*k),4,modg);
+  cv::exp(-3.315/modg, dst);
+  dst = 1.0 - dst;
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function computes a good empirical value for the k contrast factor
+ * given an input image, the percentile (0-1), the gradient scale and the number of
+ * bins in the histogram
+ * @param img Input image
+ * @param perc Percentile of the image gradient histogram (0-1)
+ * @param gscale Scale for computing the image gradient histogram
+ * @param nbins Number of histogram bins
+ * @param ksize_x Kernel size in X-direction (horizontal) for the Gaussian smoothing kernel
+ * @param ksize_y Kernel size in Y-direction (vertical) for the Gaussian smoothing kernel
+ * @return k contrast factor
+ */
+float compute_k_percentile(const cv::Mat& img, float perc, float gscale,
+                           int nbins, int ksize_x, int ksize_y) {
+
+  int nbin = 0, nelements = 0, nthreshold = 0, k = 0;
+  float kperc = 0.0, modg = 0.0, lx = 0.0, ly = 0.0;
+  float npoints = 0.0;
+  float hmax = 0.0;
+
+  // Create the array for the histogram
+  float *hist = new float[nbins];
+
+  // Create the matrices
+  Mat gaussian = Mat::zeros(img.rows,img.cols,CV_32F);
+  Mat Lx = Mat::zeros(img.rows,img.cols,CV_32F);
+  Mat Ly = Mat::zeros(img.rows,img.cols,CV_32F);
+
+  // Set the histogram to zero, just in case
+  for (int i = 0; i < nbins; i++) {
+    hist[i] = 0.0;
+  }
+
+  // Perform the Gaussian convolution
+  gaussian_2D_convolution(img,gaussian,ksize_x,ksize_y,gscale);
+
+  // Compute the Gaussian derivatives Lx and Ly
+  Scharr(gaussian,Lx,CV_32F,1,0,1,0,cv::BORDER_DEFAULT);
+  Scharr(gaussian,Ly,CV_32F,0,1,1,0,cv::BORDER_DEFAULT);
+
+  // Skip the borders for computing the histogram
+  for (int i = 1; i < gaussian.rows-1; i++) {
+    for (int j = 1; j < gaussian.cols-1; j++) {
+      lx = *(Lx.ptr<float>(i)+j);
+      ly = *(Ly.ptr<float>(i)+j);
+      modg = sqrt(lx*lx + ly*ly);
+
+      // Get the maximum
+      if (modg > hmax) {
+        hmax = modg;
+      }
+    }
+  }
+
+  // Skip the borders for computing the histogram
+  for (int i = 1; i < gaussian.rows-1; i++) {
+    for (int j = 1; j < gaussian.cols-1; j++) {
+      lx = *(Lx.ptr<float>(i)+j);
+      ly = *(Ly.ptr<float>(i)+j);
+      modg = sqrt(lx*lx + ly*ly);
+
+      // Find the correspondent bin
+      if (modg != 0.0) {
+        nbin = floor(nbins*(modg/hmax));
+
+        if (nbin == nbins) {
+          nbin--;
+        }
+
+        hist[nbin]++;
+        npoints++;
+      }
+    }
+  }
+
+  // Now find the perc of the histogram percentile
+  nthreshold = (size_t)(npoints*perc);
+
+
+  for (k = 0; nelements < nthreshold && k < nbins; k++) {
+    nelements = nelements + hist[k];
+  }
+
+  if (nelements < nthreshold)  {
+    kperc = 0.03;
+  }
+  else {
+    kperc = hmax*((float)(k)/(float)nbins);
+  }
+
+  delete hist;
+  return kperc;
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function computes Scharr image derivatives
+ * @param src Input image
+ * @param dst Output image
+ * @param xorder Derivative order in X-direction (horizontal)
+ * @param yorder Derivative order in Y-direction (vertical)
+ * @param scale Scale factor or derivative size
+ */
+void compute_scharr_derivatives(const cv::Mat& src, cv::Mat& dst,
+                                int xorder, int yorder, int scale) {
+  Mat kx, ky;
+  compute_derivative_kernels(kx,ky,xorder,yorder,scale);
+  sepFilter2D(src,dst,CV_32F,kx,ky);
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief Compute derivative kernels for sizes different than 3
+ * @param _kx Horizontal kernel values
+ * @param _ky Vertical kernel values
+ * @param dx Derivative order in X-direction (horizontal)
+ * @param dy Derivative order in Y-direction (vertical)
+ * @param scale_ Scale factor or derivative size
+ */
+void compute_derivative_kernels(cv::OutputArray _kx, cv::OutputArray _ky,
+                                int dx, int dy, int scale) {
+
+  int ksize = 3 + 2*(scale-1);
+
+  // The standard Scharr kernel
+  if (scale == 1) {
+    getDerivKernels(_kx,_ky,dx,dy,0,true,CV_32F);
+    return;
+  }
+
+  _kx.create(ksize,1,CV_32F,-1,true);
+  _ky.create(ksize,1,CV_32F,-1,true);
+  Mat kx = _kx.getMat();
+  Mat ky = _ky.getMat();
+
+  float w = 10.0/3.0;
+  float norm = 1.0/(2.0*scale*(w+2.0));
+
+  for (int k = 0; k < 2; k++) {
+    Mat* kernel = k == 0 ? &kx : &ky;
+    int order = k == 0 ? dx : dy;
+	std::vector<float> kerI(ksize);
+
+    for (int t=0; t<ksize; t++) {
+      kerI[t] = 0;
+    }
+
+    if (order == 0) {
+      kerI[0] = norm, kerI[ksize/2] = w*norm, kerI[ksize-1] = norm;
+    }
+    else if (order == 1) {
+      kerI[0] = -1, kerI[ksize/2] = 0, kerI[ksize-1] = 1;
+    }
+
+    Mat temp(kernel->rows,kernel->cols,CV_32F,&kerI[0]);
+    temp.copyTo(*kernel);
+  }
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function performs a scalar non-linear diffusion step
+ * @param Ld2 Output image in the evolution
+ * @param c Conductivity image
+ * @param Lstep Previous image in the evolution
+ * @param stepsize The step size in time units
+ * @note Forward Euler Scheme 3x3 stencil
+ * The function c is a scalar value that depends on the gradient norm
+ * dL_by_ds = d(c dL_by_dx)_by_dx + d(c dL_by_dy)_by_dy
+ */
+void nld_step_scalar(cv::Mat& Ld, const cv::Mat& c, cv::Mat& Lstep, float stepsize) {
+
+#ifdef _OPENMP
+#pragma omp parallel for schedule(dynamic)
+#endif
+  for (int i = 1; i < Lstep.rows-1; i++) {
+    for (int j = 1; j < Lstep.cols-1; j++) {
+      float xpos = ((*(c.ptr<float>(i)+j))+(*(c.ptr<float>(i)+j+1)))*((*(Ld.ptr<float>(i)+j+1))-(*(Ld.ptr<float>(i)+j)));
+      float xneg = ((*(c.ptr<float>(i)+j-1))+(*(c.ptr<float>(i)+j)))*((*(Ld.ptr<float>(i)+j))-(*(Ld.ptr<float>(i)+j-1)));
+      float ypos = ((*(c.ptr<float>(i)+j))+(*(c.ptr<float>(i+1)+j)))*((*(Ld.ptr<float>(i+1)+j))-(*(Ld.ptr<float>(i)+j)));
+      float yneg = ((*(c.ptr<float>(i-1)+j))+(*(c.ptr<float>(i)+j)))*((*(Ld.ptr<float>(i)+j))-(*(Ld.ptr<float>(i-1)+j)));
+      *(Lstep.ptr<float>(i)+j) = 0.5*stepsize*(xpos-xneg + ypos-yneg);
+    }
+  }
+
+  for (int j = 1; j < Lstep.cols-1; j++) {
+    float xpos = ((*(c.ptr<float>(0)+j))+(*(c.ptr<float>(0)+j+1)))*((*(Ld.ptr<float>(0)+j+1))-(*(Ld.ptr<float>(0)+j)));
+    float xneg = ((*(c.ptr<float>(0)+j-1))+(*(c.ptr<float>(0)+j)))*((*(Ld.ptr<float>(0)+j))-(*(Ld.ptr<float>(0)+j-1)));
+    float ypos = ((*(c.ptr<float>(0)+j))+(*(c.ptr<float>(1)+j)))*((*(Ld.ptr<float>(1)+j))-(*(Ld.ptr<float>(0)+j)));
+    float yneg = ((*(c.ptr<float>(0)+j))+(*(c.ptr<float>(0)+j)))*((*(Ld.ptr<float>(0)+j))-(*(Ld.ptr<float>(0)+j)));
+    *(Lstep.ptr<float>(0)+j) = 0.5*stepsize*(xpos-xneg + ypos-yneg);
+  }
+
+  for (int j = 1; j < Lstep.cols-1; j++) {
+    float xpos = ((*(c.ptr<float>(Lstep.rows-1)+j))+(*(c.ptr<float>(Lstep.rows-1)+j+1)))*((*(Ld.ptr<float>(Lstep.rows-1)+j+1))-(*(Ld.ptr<float>(Lstep.rows-1)+j)));
+    float xneg = ((*(c.ptr<float>(Lstep.rows-1)+j-1))+(*(c.ptr<float>(Lstep.rows-1)+j)))*((*(Ld.ptr<float>(Lstep.rows-1)+j))-(*(Ld.ptr<float>(Lstep.rows-1)+j-1)));
+    float ypos = ((*(c.ptr<float>(Lstep.rows-1)+j))+(*(c.ptr<float>(Lstep.rows-1)+j)))*((*(Ld.ptr<float>(Lstep.rows-1)+j))-(*(Ld.ptr<float>(Lstep.rows-1)+j)));
+    float yneg = ((*(c.ptr<float>(Lstep.rows-2)+j))+(*(c.ptr<float>(Lstep.rows-1)+j)))*((*(Ld.ptr<float>(Lstep.rows-1)+j))-(*(Ld.ptr<float>(Lstep.rows-2)+j)));
+    *(Lstep.ptr<float>(Lstep.rows-1)+j) = 0.5*stepsize*(xpos-xneg + ypos-yneg);
+  }
+
+  for (int i = 1; i < Lstep.rows-1; i++) {
+    float xpos = ((*(c.ptr<float>(i)))+(*(c.ptr<float>(i)+1)))*((*(Ld.ptr<float>(i)+1))-(*(Ld.ptr<float>(i))));
+    float xneg = ((*(c.ptr<float>(i)))+(*(c.ptr<float>(i))))*((*(Ld.ptr<float>(i)))-(*(Ld.ptr<float>(i))));
+    float ypos = ((*(c.ptr<float>(i)))+(*(c.ptr<float>(i+1))))*((*(Ld.ptr<float>(i+1)))-(*(Ld.ptr<float>(i))));
+    float yneg = ((*(c.ptr<float>(i-1)))+(*(c.ptr<float>(i))))*((*(Ld.ptr<float>(i)))-(*(Ld.ptr<float>(i-1))));
+    *(Lstep.ptr<float>(i)) = 0.5*stepsize*(xpos-xneg + ypos-yneg);
+  }
+
+  for (int i = 1; i < Lstep.rows-1; i++) {
+    float xpos = ((*(c.ptr<float>(i)+Lstep.cols-1))+(*(c.ptr<float>(i)+Lstep.cols-1)))*((*(Ld.ptr<float>(i)+Lstep.cols-1))-(*(Ld.ptr<float>(i)+Lstep.cols-1)));
+    float xneg = ((*(c.ptr<float>(i)+Lstep.cols-2))+(*(c.ptr<float>(i)+Lstep.cols-1)))*((*(Ld.ptr<float>(i)+Lstep.cols-1))-(*(Ld.ptr<float>(i)+Lstep.cols-2)));
+    float ypos = ((*(c.ptr<float>(i)+Lstep.cols-1))+(*(c.ptr<float>(i+1)+Lstep.cols-1)))*((*(Ld.ptr<float>(i+1)+Lstep.cols-1))-(*(Ld.ptr<float>(i)+Lstep.cols-1)));
+    float yneg = ((*(c.ptr<float>(i-1)+Lstep.cols-1))+(*(c.ptr<float>(i)+Lstep.cols-1)))*((*(Ld.ptr<float>(i)+Lstep.cols-1))-(*(Ld.ptr<float>(i-1)+Lstep.cols-1)));
+    *(Lstep.ptr<float>(i)+Lstep.cols-1) = 0.5*stepsize*(xpos-xneg + ypos-yneg);
+  }
+
+  Ld = Ld + Lstep;
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function checks if a given pixel is a maximum in a local neighbourhood
+ * @param img Input image where we will perform the maximum search
+ * @param dsize Half size of the neighbourhood
+ * @param value Response value at (x,y) position
+ * @param row Image row coordinate
+ * @param col Image column coordinate
+ * @param same_img Flag to indicate if the image value at (x,y) is in the input image
+ * @return 1->is maximum, 0->otherwise
+ */
+bool check_maximum_neighbourhood(const cv::Mat& img, int dsize, float value,
+                                 int row, int col, bool same_img) {
+
+  bool response = true;
+
+  for (int i = row-dsize; i <= row+dsize; i++) {
+    for (int j = col-dsize; j <= col+dsize; j++) {
+      if (i >= 0 && i < img.rows && j >= 0 && j < img.cols) {
+        if (same_img == true) {
+          if (i != row || j != col) {
+            if ((*(img.ptr<float>(i)+j)) > value) {
+              response = false;
+              return response;
+            }
+          }
+        }
+        else {
+          if ((*(img.ptr<float>(i)+j)) > value) {
+            response = false;
+            return response;
+          }
+        }
+      }
+    }
+  }
+
+  return response;
+}

modules/features2d/src/kaze/nldiffusion_functions.h

@@ -0,0 +1,51 @@
+
+/**
+ * @file nldiffusion_functions.h
+ * @brief Functions for non-linear diffusion applications:
+ * 2D Gaussian Derivatives
+ * Perona and Malik conductivity equations
+ * Perona and Malik evolution
+ * @date Dec 27, 2011
+ * @author Pablo F. Alcantarilla
+ */
+
+#ifndef NLDIFFUSION_FUNCTIONS_H_
+#define NLDIFFUSION_FUNCTIONS_H_
+
+//******************************************************************************
+//******************************************************************************
+
+// Includes
+#include "config.h"
+
+//*************************************************************************************
+//*************************************************************************************
+
+// Gaussian 2D convolution
+void gaussian_2D_convolution(const cv::Mat& src, cv::Mat& dst,
+                             int ksize_x, int ksize_y, float sigma);
+
+// Diffusivity functions
+void pm_g1(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k);
+void pm_g2(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k);
+void weickert_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k);
+float compute_k_percentile(const cv::Mat& img, float perc, float gscale,
+                           int nbins, int ksize_x, int ksize_y);
+
+// Image derivatives
+void compute_scharr_derivatives(const cv::Mat& src, cv::Mat& dst,
+                                int xorder, int yorder, int scale);
+void compute_derivative_kernels(cv::OutputArray _kx, cv::OutputArray _ky,
+                                int dx, int dy, int scale);
+
+// Nonlinear diffusion filtering scalar step
+void nld_step_scalar(cv::Mat& Ld, const cv::Mat& c, cv::Mat& Lstep, float stepsize);
+
+// For non-maxima suppresion
+bool check_maximum_neighbourhood(const cv::Mat& img, int dsize, float value,
+                                 int row, int col, bool same_img);
+
+//*************************************************************************************
+//*************************************************************************************
+
+#endif // NLDIFFUSION_FUNCTIONS_H_

modules/features2d/src/kaze/utils.cpp

@@ -0,0 +1,92 @@
+
+//=============================================================================
+//
+// utils.cpp
+// Author: Pablo F. Alcantarilla
+// Institution: University d'Auvergne
+// Address: Clermont Ferrand, France
+// Date: 29/12/2011
+// Email: pablofdezalc@gmail.com
+//
+// KAZE Features Copyright 2012, Pablo F. Alcantarilla
+// All Rights Reserved
+// See LICENSE for the license information
+//=============================================================================
+
+/**
+ * @file utils.cpp
+ * @brief Some useful functions
+ * @date Dec 29, 2011
+ * @author Pablo F. Alcantarilla
+ */
+
+#include "utils.h"
+
+using namespace std;
+using namespace cv;
+
+//*************************************************************************************
+//*************************************************************************************
+
+/**
+ * @brief This function copies the input image and converts the scale of the copied
+ * image prior visualization
+ * @param src Input image
+ * @param dst Output image
+ */
+void copy_and_convert_scale(const cv::Mat& src, cv::Mat& dst) {
+
+  float min_val = 0, max_val = 0;
+
+  src.copyTo(dst);
+  compute_min_32F(dst,min_val);
+
+  dst = dst - min_val;
+
+  compute_max_32F(dst,max_val);
+  dst = dst / max_val;
+}
+
+//*************************************************************************************
+//*************************************************************************************
+
+/*
+void show_input_options_help(int example) {
+
+  fflush(stdout);
+
+  cout << endl;
+  cout << endl;
+  cout << "KAZE Features" << endl;
+  cout << "***********************************************************" << endl;
+  cout << "For running the program you need to type in the command line the following arguments: " << endl;
+
+  if (example == 0) {
+    cout << "./kaze_features img.jpg [options]" << endl;
+  }
+  else if (example == 1) {
+    cout << "./kaze_match img1.jpg img2.pgm homography.txt [options]" << endl;
+  }
+  else if (example == 2) {
+    cout << "./kaze_compare img1.jpg img2.pgm homography.txt [options]" << endl;
+  }
+
+  cout << endl;
+  cout << "The options are not mandatory. In case you do not specify additional options, default arguments will be used" << endl << endl;
+  cout << "Here is a description of the additional options: " << endl;
+  cout << "--verbose " << "\t\t if verbosity is required" << endl;
+  cout << "--help" << "\t\t for showing the command line options" << endl;
+  cout << "--soffset" << "\t\t the base scale offset (sigma units)" << endl;
+  cout << "--omax" << "\t\t maximum octave evolution of the image 2^sigma (coarsest scale)" << endl;
+  cout << "--nsublevels" << "\t\t number of sublevels per octave" << endl;
+  cout << "--dthreshold" << "\t\t Feature detector threshold response for accepting points (0.001 can be a good value)" << endl;
+  cout << "--descriptor" << "\t\t Descriptor Type 0 -> SURF, 1 -> M-SURF, 2 -> G-SURF" << endl;
+  cout << "--use_fed" "\t\t 1 -> Use FED, 0 -> Use AOS for the nonlinear diffusion filtering" << endl;
+  cout << "--upright" << "\t\t 0 -> Rotation Invariant, 1 -> No Rotation Invariant" << endl;
+  cout << "--extended" << "\t\t 0 -> Normal Descriptor (64), 1 -> Extended Descriptor (128)" << endl;
+  cout << "--output keypoints.txt" << "\t\t For saving the detected keypoints into a .txt file" << endl;
+  cout << "--save_scale_space" << "\t\t 1 in case we want to save the nonlinear scale space images. 0 otherwise" << endl;
+  cout << "--show_results" << "\t\t 1 in case we want to show detection results. 0 otherwise" << endl;
+  cout << endl;
+}
+*/

modules/features2d/src/kaze/utils.h

@@ -0,0 +1,41 @@
+
+/**
+ * @file utils.h
+ * @brief Some useful functions
+ * @date Dec 29, 2011
+ * @author Pablo F. Alcantarilla
+ */
+
+#ifndef UTILS_H_
+#define UTILS_H_
+
+//******************************************************************************
+//******************************************************************************
+
+// OPENCV Includes
+#include "precomp.hpp"
+
+// System Includes
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <cstdlib>
+#include <string>
+#include <vector>
+#include <fstream>
+#include <assert.h>
+#include <math.h>
+
+//*************************************************************************************
+//*************************************************************************************
+
+// Declaration of Functions
+void compute_min_32F(const cv::Mat& src, float& value);
+void compute_max_32F(const cv::Mat& src, float& value);
+void convert_scale(cv::Mat& src);
+void copy_and_convert_scale(const cv::Mat &src, cv::Mat& dst);
+
+//*************************************************************************************
+//*************************************************************************************
+
+#endif // UTILS_H_