Merge pull request #8951 from hrnr:akaze_part2

[GSOC] Speeding-up AKAZE, part #2 (#8951)

* feature2d: instrument more functions used in AKAZE

* rework Compute_Determinant_Hessian_Response

* this takes 84% of time of Feature_Detection
* run everything in parallel
* compute Scharr kernels just once
* compute sigma more efficiently
* allocate all matrices in evolution without zeroing

* features2d: add one bigger image to tests

* now test have images: 600x768, 900x600 and 1385x700 to cover different resolutions

* explicitly zero Lx and Ly

* add Lflow and Lstep to evolution as in original AKAZE code

* reworked computing keypoints orientation

integrated faster function from https://github.com/h2suzuki/fast_akaze

* use standard fastAtan2 instead of getAngle

* compute keypoints orientation in parallel

* fix visual studio warnings

* replace some wrapped functions with direct calls to OpenCV functions

* improved readability for people familiar with opencv
* do not same image twice in base level

* rework diffusity stencil

* use one pass stencil for diffusity from https://github.com/h2suzuki/fast_akaze
* improve locality in Create_Scale_Space

* always compute determinat od hessian and spacial derivatives

* this needs to be computed always as we need derivatives while computing descriptors
* fixed tests of AKAZE with KAZE descriptors which have been affected by this

Currently it computes all first and second order derivatives together and the determiant of the hessian. For descriptors it would be enough to compute just first order derivates, but it is not probably worth it optimize for scenario where descriptors and keypoints are computed separately, since it is already very inefficient. When computing keypoint and descriptors together it is faster to do it the current way (preserves locality).

* parallelize non linear diffusion computation

* do multiplication right in the nlp diffusity kernel

* rework kfactor computation

* get rid of sharing buffers when creating scale space pyramid, the performace impact is neglegible

* features2d: initialize TBB scheduler in perf tests

* ensures more stable output
* more reasonable profiles, since the first call of parallel_for_ is not getting big performace hit

* compute_kfactor: interleave finding of maximum and computing distance

* no need to go twice through the data

* start to use UMats in AKAZE to leverage OpenCl in the future

* fixed bug that prevented computing determinant for scale pyramid of size 1 (just the base image)
* all descriptors now support writing to uninitialized memory
* use InputArray and OutputArray for input image and descriptors, allows to make use UMAt that user passes to us

* enable use of all existing ocl paths in AKAZE

* all parts that uses ocl-enabled functions should use ocl by now

* imgproc: fix dispatching of IPP version when OCL is disabled

* when OCL is disabled IPP version should be always prefered (even when the dst is UMat)

* get rid of copy in DeterminantHessian response

* this slows CPU version considerably
* do no run in parallel when running with OCL

* store derivations as UMat in pyramid

* enables OCL path computing of determint hessian
* will allow to compute descriptors on GPU in the future

* port diffusivity to OCL

* diffusivity itself is not a blocker, but this saves us downloading and uploading derivations

* implement kernel for nonlinear scalar diffusion step

* download the pyramid from GPU just once

we don't want to downlaod matrices ad hoc from gpu when the function in AKAZE needs it. There is a HUGE mapping overhead and without shared memory support a LOT of unnecessary transfers.

This maps/downloads matrices just once.

* fix bug with uninitialized values in non linear diffusion

* this was causing spurious segfaults in stitching tests due to propagation of NaNs
* added new test, which checks for NaNs (added new debug asserts for NaNs)
* valgrind now says everything is ok

* add nonlinear diffusion step OCL implementation

* Lt in pyramid changed to UMat, it will be downlaoded from GPU along with Lx, Ly
* fix bug in pm_g2 kernel. OpenCV mangles dimensions passed to OpenCL, so we need to check for boundaries in each OCL kernel.

* port computing of determinant to OCL

* computing of determinant is not a blocker, but with this change we don't need to download all spatial derivatives to CPU, we only download determinant
* make Ldet in the pyramid UMat, download it from CPU together with the other parts of the pyramid
* add profiling macros

* fix visual studio warning

* instrument non_linear_diffusion

* remove changes I have made to TEvolution

* TEvolution is used only in KAZE now

* Revert "features2d: initialize TBB scheduler in perf tests"

This reverts commit ba81e2a711ae009ce3c5459775627b6423112669.

modules/features2d/perf/perf_feature2d.hpp

@@ -35,7 +35,8 @@ typedef perf::TestBaseWithParam<Feature2DType_String_t> feature2d;
 
 #define TEST_IMAGES testing::Values(\
     "cv/detectors_descriptors_evaluation/images_datasets/leuven/img1.png",\
-    "stitching/a3.png")
+    "stitching/a3.png", \
+    "stitching/s2.jpg")
 
 static inline Ptr<Feature2D> getFeature2D(Feature2DType type)
 {

modules/features2d/src/akaze.cpp

@@ -210,13 +210,10 @@ namespace cv
 
             if( descriptors.needed() )
             {
-                Mat desc;
-                impl.Compute_Descriptors(keypoints, desc);
-                // TODO optimize this copy
-                desc.copyTo(descriptors);
+                impl.Compute_Descriptors(keypoints, descriptors);
 
-                CV_Assert((!desc.rows || desc.cols == descriptorSize()));
-                CV_Assert((!desc.rows || (desc.type() == descriptorType())));
+                CV_Assert((descriptors.empty() || descriptors.cols() == descriptorSize()));
+                CV_Assert((descriptors.empty() || (descriptors.type() == descriptorType())));
             }
         }
 

modules/features2d/src/kaze/AKAZEFeatures.cpp

@@ -11,6 +11,7 @@
 #include "fed.h"
 #include "nldiffusion_functions.h"
 #include "utils.h"
+#include "opencl_kernels_features2d.hpp"
 
 #include <iostream>
 
@@ -43,6 +44,7 @@ AKAZEFeatures::AKAZEFeatures(const AKAZEOptions& options) : options_(options) {
  * @brief This method allocates the memory for the nonlinear diffusion evolution
  */
 void AKAZEFeatures::Allocate_Memory_Evolution(void) {
+  CV_INSTRUMENT_REGION()
 
   float rfactor = 0.0f;
   int level_height = 0, level_width = 0;
@@ -60,20 +62,14 @@ void AKAZEFeatures::Allocate_Memory_Evolution(void) {
     }
 
     for (int j = 0; j < options_.nsublevels; j++) {
-      TEvolution step;
-      step.Lx = Mat::zeros(level_height, level_width, CV_32F);
-      step.Ly = Mat::zeros(level_height, level_width, CV_32F);
-      step.Lxx = Mat::zeros(level_height, level_width, CV_32F);
-      step.Lxy = Mat::zeros(level_height, level_width, CV_32F);
-      step.Lyy = Mat::zeros(level_height, level_width, CV_32F);
-      step.Lt = Mat::zeros(level_height, level_width, CV_32F);
-      step.Ldet = Mat::zeros(level_height, level_width, CV_32F);
-      step.Lsmooth = Mat::zeros(level_height, level_width, CV_32F);
+      Evolution step;
+      step.size = Size(level_width, level_height);
       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.sigma_size = fRound(step.esigma * options_.derivative_factor / power);  // In fact sigma_size only depends on j
+      step.etime = 0.5f * (step.esigma * step.esigma);
       step.octave = i;
       step.sublevel = j;
+      step.octave_ratio = (float)power;
       evolution_.push_back(step);
     }
   }
@@ -92,74 +88,393 @@ void AKAZEFeatures::Allocate_Memory_Evolution(void) {
 }
 
 /* ************************************************************************* */
+/**
+ * @brief Computes kernel size for Gaussian smoothing if the image
+ * @param sigma Kernel standard deviation
+ * @returns kernel size
+ */
+static inline int getGaussianKernelSize(float sigma) {
+  // Compute an appropriate kernel size according to the specified sigma
+  int ksize = (int)ceil(2.0f*(1.0f + (sigma - 0.8f) / (0.3f)));
+  ksize |= 1; // kernel should be odd
+  return ksize;
+}
+
+/* ************************************************************************* */
+/**
+* @brief This function computes a scalar non-linear diffusion step
+* @param Lt Base image in the evolution
+* @param Lf Conductivity image
+* @param Lstep Output image that gives the difference between the current
+* Ld and the next Ld being evolved
+* @param row_begin row where to start
+* @param row_end last row to fill exclusive. the range is [row_begin, row_end).
+* @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
+*/
+static inline void
+nld_step_scalar_one_lane(const Mat& Lt, const Mat& Lf, Mat& Lstep, float step_size, int row_begin, int row_end)
+{
+  CV_INSTRUMENT_REGION()
+  /* The labeling scheme for this five star stencil:
+   [    a    ]
+   [ -1 c +1 ]
+   [    b    ]
+   */
+
+  Lstep.create(Lt.size(), Lt.type());
+  const int cols = Lt.cols - 2;
+  int row = row_begin;
+
+  const float *lt_a, *lt_c, *lt_b;
+  const float *lf_a, *lf_c, *lf_b;
+  float *dst;
+  float step_r = 0.f;
+
+  // Process the top row
+  if (row == 0) {
+    lt_c = Lt.ptr<float>(0) + 1;  /* Skip the left-most column by +1 */
+    lf_c = Lf.ptr<float>(0) + 1;
+    lt_b = Lt.ptr<float>(1) + 1;
+    lf_b = Lf.ptr<float>(1) + 1;
+
+    // fill the corner to prevent uninitialized values
+    dst = Lstep.ptr<float>(0);
+    dst[0] = 0.0f;
+    ++dst;
+
+    for (int j = 0; j < cols; j++) {
+      step_r = (lf_c[j] + lf_c[j + 1])*(lt_c[j + 1] - lt_c[j]) +
+               (lf_c[j] + lf_c[j - 1])*(lt_c[j - 1] - lt_c[j]) +
+               (lf_c[j] + lf_b[j    ])*(lt_b[j    ] - lt_c[j]);
+      dst[j] = step_r * step_size;
+    }
+
+    // fill the corner to prevent uninitialized values
+    dst[cols] = 0.0f;
+    ++row;
+  }
+
+  // Process the middle rows
+  int middle_end = std::min(Lt.rows - 1, row_end);
+  for (; row < middle_end; ++row)
+  {
+    lt_a = Lt.ptr<float>(row - 1);
+    lf_a = Lf.ptr<float>(row - 1);
+    lt_c = Lt.ptr<float>(row    );
+    lf_c = Lf.ptr<float>(row    );
+    lt_b = Lt.ptr<float>(row + 1);
+    lf_b = Lf.ptr<float>(row + 1);
+    dst = Lstep.ptr<float>(row);
+
+    // The left-most column
+    step_r = (lf_c[0] + lf_c[1])*(lt_c[1] - lt_c[0]) +
+             (lf_c[0] + lf_b[0])*(lt_b[0] - lt_c[0]) +
+             (lf_c[0] + lf_a[0])*(lt_a[0] - lt_c[0]);
+    dst[0] = step_r * step_size;
+
+    lt_a++; lt_c++; lt_b++;
+    lf_a++; lf_c++; lf_b++;
+    dst++;
+
+    // The middle columns
+    for (int j = 0; j < cols; j++)
+    {
+      step_r = (lf_c[j] + lf_c[j + 1])*(lt_c[j + 1] - lt_c[j]) +
+               (lf_c[j] + lf_c[j - 1])*(lt_c[j - 1] - lt_c[j]) +
+               (lf_c[j] + lf_b[j    ])*(lt_b[j    ] - lt_c[j]) +
+               (lf_c[j] + lf_a[j    ])*(lt_a[j    ] - lt_c[j]);
+      dst[j] = step_r * step_size;
+    }
+
+    // The right-most column
+    step_r = (lf_c[cols] + lf_c[cols - 1])*(lt_c[cols - 1] - lt_c[cols]) +
+             (lf_c[cols] + lf_b[cols    ])*(lt_b[cols    ] - lt_c[cols]) +
+             (lf_c[cols] + lf_a[cols    ])*(lt_a[cols    ] - lt_c[cols]);
+    dst[cols] = step_r * step_size;
+  }
+
+  // Process the bottom row (row == Lt.rows - 1)
+  if (row_end == Lt.rows) {
+    lt_a = Lt.ptr<float>(row - 1) + 1;  /* Skip the left-most column by +1 */
+    lf_a = Lf.ptr<float>(row - 1) + 1;
+    lt_c = Lt.ptr<float>(row    ) + 1;
+    lf_c = Lf.ptr<float>(row    ) + 1;
+
+    // fill the corner to prevent uninitialized values
+    dst = Lstep.ptr<float>(row);
+    dst[0] = 0.0f;
+    ++dst;
+
+    for (int j = 0; j < cols; j++) {
+      step_r = (lf_c[j] + lf_c[j + 1])*(lt_c[j + 1] - lt_c[j]) +
+               (lf_c[j] + lf_c[j - 1])*(lt_c[j - 1] - lt_c[j]) +
+               (lf_c[j] + lf_a[j    ])*(lt_a[j    ] - lt_c[j]);
+      dst[j] = step_r * step_size;
+    }
+
+    // fill the corner to prevent uninitialized values
+    dst[cols] = 0.0f;
+  }
+}
+
+class NonLinearScalarDiffusionStep : public ParallelLoopBody
+{
+public:
+  NonLinearScalarDiffusionStep(const Mat& Lt, const Mat& Lf, Mat& Lstep, float step_size)
+    : Lt_(&Lt), Lf_(&Lf), Lstep_(&Lstep), step_size_(step_size)
+  {}
+
+  void operator()(const Range& range) const
+  {
+    nld_step_scalar_one_lane(*Lt_, *Lf_, *Lstep_, step_size_, range.start, range.end);
+  }
+
+private:
+  const Mat* Lt_;
+  const Mat* Lf_;
+  Mat* Lstep_;
+  float step_size_;
+};
+
+#ifdef HAVE_OPENCL
+static inline bool
+ocl_non_linear_diffusion_step(const UMat& Lt, const UMat& Lf, UMat& Lstep, float step_size)
+{
+  if(!Lt.isContinuous())
+    return false;
+
+  size_t globalSize[] = {(size_t)Lt.cols, (size_t)Lt.rows};
+
+  ocl::Kernel ker("AKAZE_nld_step_scalar", ocl::features2d::akaze_oclsrc);
+  if( ker.empty() )
+    return false;
+
+  return ker.args(
+    ocl::KernelArg::ReadOnly(Lt),
+    ocl::KernelArg::PtrReadOnly(Lf),
+    ocl::KernelArg::PtrWriteOnly(Lstep),
+    step_size).run(2, globalSize, 0, true);
+}
+#endif // HAVE_OPENCL
+
+static inline void
+non_linear_diffusion_step(const UMat& Lt, const UMat& Lf, UMat& Lstep, float step_size)
+{
+  CV_INSTRUMENT_REGION()
+
+  Lstep.create(Lt.size(), Lt.type());
+
+  CV_OCL_RUN(true, ocl_non_linear_diffusion_step(Lt, Lf, Lstep, step_size));
+
+  // when on CPU UMats should be already allocated on CPU so getMat here is basicallly no-op
+  Mat Mstep = Lstep.getMat(ACCESS_WRITE);
+  parallel_for_(Range(0, Lt.rows), NonLinearScalarDiffusionStep(Lt.getMat(ACCESS_READ),
+    Lf.getMat(ACCESS_READ), Mstep, step_size));
+}
+
+/**
+ * @brief This function computes a good empirical value for the k contrast factor
+ * given two gradient images, the percentile (0-1), the temporal storage to hold
+ * gradient norms and the histogram bins
+ * @param Lx Horizontal gradient of the input image
+ * @param Ly Vertical gradient of the input image
+ * @param nbins Number of histogram bins
+ * @return k contrast factor
+ */
+static inline float
+compute_kcontrast(const cv::Mat& Lx, const cv::Mat& Ly, float perc, int nbins)
+{
+  CV_INSTRUMENT_REGION()
+
+  CV_Assert(nbins > 2);
+  CV_Assert(!Lx.empty());
+
+  // temporary square roots of dot product
+  Mat modgs (Lx.rows - 2, Lx.cols - 2, CV_32F);
+  const int total = modgs.cols * modgs.rows;
+  float *modg = modgs.ptr<float>();
+  float hmax = 0.0f;
+
+  for (int i = 1; i < Lx.rows - 1; i++) {
+    const float *lx = Lx.ptr<float>(i) + 1;
+    const float *ly = Ly.ptr<float>(i) + 1;
+    const int cols = Lx.cols - 2;
+
+    for (int j = 0; j < cols; j++) {
+      float dist = sqrtf(lx[j] * lx[j] + ly[j] * ly[j]);
+      *modg++ = dist;
+      hmax = std::max(hmax, dist);
+    }
+  }
+  modg = modgs.ptr<float>();
+
+  if (hmax == 0.0f)
+    return 0.03f;  // e.g. a blank image
+
+  // Compute the bin numbers: the value range [0, hmax] -> [0, nbins-1]
+  modgs *= (nbins - 1) / hmax;
+
+  // Count up histogram
+  std::vector<int> hist(nbins, 0);
+  for (int i = 0; i < total; i++)
+    hist[(int)modg[i]]++;
+
+  // Now find the perc of the histogram percentile
+  const int nthreshold = (int)((total - hist[0]) * perc);  // Exclude hist[0] as background
+  int nelements = 0;
+  for (int k = 1; k < nbins; k++) {
+    if (nelements >= nthreshold)
+        return (float)hmax * k / nbins;
+
+    nelements += hist[k];
+  }
+
+  return 0.03f;
+}
+
+#ifdef HAVE_OPENCL
+static inline bool
+ocl_pm_g2(const UMat& Lx, const UMat& Ly, UMat& Lflow, float kcontrast)
+{
+  int total = Lx.rows * Lx.cols;
+  size_t globalSize[] = {(size_t)total};
+
+  ocl::Kernel ker("AKAZE_pm_g2", ocl::features2d::akaze_oclsrc);
+  if( ker.empty() )
+    return false;
+
+  return ker.args(
+    ocl::KernelArg::PtrReadOnly(Lx),
+    ocl::KernelArg::PtrReadOnly(Ly),
+    ocl::KernelArg::PtrWriteOnly(Lflow),
+    kcontrast, total).run(1, globalSize, 0, true);
+}
+#endif // HAVE_OPENCL
+
+static inline void
+compute_diffusivity(const UMat& Lx, const UMat& Ly, UMat& Lflow, float kcontrast, int diffusivity)
+{
+  CV_INSTRUMENT_REGION()
+
+  Lflow.create(Lx.size(), Lx.type());
+
+  switch (diffusivity) {
+    case KAZE::DIFF_PM_G1:
+      pm_g1(Lx, Ly, Lflow, kcontrast);
+    break;
+    case KAZE::DIFF_PM_G2:
+      CV_OCL_RUN(true, ocl_pm_g2(Lx, Ly, Lflow, kcontrast));
+      pm_g2(Lx, Ly, Lflow, kcontrast);
+    break;
+    case KAZE::DIFF_WEICKERT:
+      weickert_diffusivity(Lx, Ly, Lflow, kcontrast);
+    break;
+    case KAZE::DIFF_CHARBONNIER:
+      charbonnier_diffusivity(Lx, Ly, Lflow, kcontrast);
+    break;
+    default:
+      CV_Error(diffusivity, "Diffusivity is not supported");
+    break;
+  }
+}
+
+/**
+ * @brief Fetches pyramid from the gpu.
+ * @details Setups mapping for matrices that might be probably on the GPU, if the
+ * code executes with OpenCL. This will setup MLx, MLy, Mdet members in the pyramid with
+ * mapping to respective UMats. This must be called before CPU-only parts of AKAZE, that work
+ * only on these Mats.
+ *
+ * This prevents mapping/unmapping overhead (and possible uploads/downloads) that would occur, if
+ * we just create Mats from UMats each time we need it later. This has devastating effects on OCL
+ * performace.
+ *
+ * @param evolution Pyramid to download
+ */
+static inline void downloadPyramid(std::vector<Evolution>& evolution)
+{
+  CV_INSTRUMENT_REGION()
+
+  for (size_t i = 0; i < evolution.size(); ++i) {
+    Evolution& e = evolution[i];
+    e.Mx = e.Lx.getMat(ACCESS_READ);
+    e.My = e.Ly.getMat(ACCESS_READ);
+    e.Mt = e.Lt.getMat(ACCESS_READ);
+    e.Mdet = e.Ldet.getMat(ACCESS_READ);
+  }
+}
+
 /**
  * @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 Mat& img)
+void AKAZEFeatures::Create_Nonlinear_Scale_Space(InputArray img)
 {
+  CV_INSTRUMENT_REGION()
   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);
+  // create first level of the evolution
+  int ksize = getGaussianKernelSize(options_.soffset);
+  GaussianBlur(img, evolution_[0].Lsmooth, Size(ksize, ksize), options_.soffset, options_.soffset, BORDER_REPLICATE);
+  evolution_[0].Lsmooth.copyTo(evolution_[0].Lt);
 
-  // Allocate memory for the flow and step images
-  Mat Lflow = Mat::zeros(evolution_[0].Lt.rows, evolution_[0].Lt.cols, CV_32F);
-  Mat Lstep = Mat::zeros(evolution_[0].Lt.rows, evolution_[0].Lt.cols, CV_32F);
+  if (evolution_.size() == 1) {
+    // we don't need to compute kcontrast factor
+    Compute_Determinant_Hessian_Response();
+    downloadPyramid(evolution_);
+    return;
+  }
+
+  // derivatives, flow and diffusion step
+  UMat Lx, Ly, Lsmooth, Lflow, Lstep;
 
-  // First compute the kcontrast factor
-  options_.kcontrast = compute_k_percentile(img, options_.kcontrast_percentile, 1.0f, options_.kcontrast_nbins, 0, 0);
+  // compute derivatives for computing k contrast
+  GaussianBlur(img, Lsmooth, Size(5, 5), 1.0f, 1.0f, BORDER_REPLICATE);
+  Scharr(Lsmooth, Lx, CV_32F, 1, 0, 1, 0, BORDER_DEFAULT);
+  Scharr(Lsmooth, Ly, CV_32F, 0, 1, 1, 0, BORDER_DEFAULT);
+  Lsmooth.release();
+  // compute the kcontrast factor
+  float kcontrast = compute_kcontrast(Lx.getMat(ACCESS_READ), Ly.getMat(ACCESS_READ),
+    options_.kcontrast_percentile, options_.kcontrast_nbins);
 
   // Now generate the rest of evolution levels
   for (size_t i = 1; i < evolution_.size(); i++) {
+    Evolution &e = evolution_[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 = Mat::zeros(evolution_[i].Lt.rows, evolution_[i].Lt.cols, CV_32F);
-      Lstep = Mat::zeros(evolution_[i].Lt.rows, evolution_[i].Lt.cols, CV_32F);
+    if (e.octave > evolution_[i - 1].octave) {
+      // new octave will be half the size
+      resize(evolution_[i - 1].Lt, e.Lt, e.size, 0, 0, INTER_AREA);
+      kcontrast *= 0.75f;
     }
     else {
-      evolution_[i - 1].Lt.copyTo(evolution_[i].Lt);
+      evolution_[i - 1].Lt.copyTo(e.Lt);
     }
 
-    gaussian_2D_convolution(evolution_[i].Lt, evolution_[i].Lsmooth, 0, 0, 1.0f);
+    GaussianBlur(e.Lt, e.Lsmooth, Size(5, 5), 1.0f, 1.0f, BORDER_REPLICATE);
 
     // 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);
+    Scharr(e.Lsmooth, Lx, CV_32F, 1, 0, 1.0, 0, BORDER_DEFAULT);
+    Scharr(e.Lsmooth, Ly, CV_32F, 0, 1, 1.0, 0, BORDER_DEFAULT);
 
     // Compute the conductivity equation
-    switch (options_.diffusivity) {
-      case KAZE::DIFF_PM_G1:
-        pm_g1(evolution_[i].Lx, evolution_[i].Ly, Lflow, options_.kcontrast);
-      break;
-      case KAZE::DIFF_PM_G2:
-        pm_g2(evolution_[i].Lx, evolution_[i].Ly, Lflow, options_.kcontrast);
-      break;
-      case KAZE::DIFF_WEICKERT:
-        weickert_diffusivity(evolution_[i].Lx, evolution_[i].Ly, Lflow, options_.kcontrast);
-      break;
-      case KAZE::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++) {
-      nld_step_scalar(evolution_[i].Lt, Lflow, Lstep, tsteps_[i - 1][j]);
+    compute_diffusivity(Lx, Ly, Lflow, kcontrast, options_.diffusivity);
+
+    // Perform Fast Explicit Diffusion on Lt
+    std::vector<float> &tsteps = tsteps_[i - 1];
+    for (size_t j = 0; j < tsteps.size(); j++) {
+      const float step_size = tsteps[j] * 0.5f;
+      non_linear_diffusion_step(e.Lt, Lflow, Lstep, step_size);
+      add(e.Lt, Lstep, e.Lt);
     }
   }
 
-  return 0;
+  Compute_Determinant_Hessian_Response();
+  downloadPyramid(evolution_);
+
+  return;
 }
 
 /* ************************************************************************* */
@@ -169,82 +484,126 @@ int AKAZEFeatures::Create_Nonlinear_Scale_Space(const Mat& img)
  */
 void AKAZEFeatures::Feature_Detection(std::vector<KeyPoint>& kpts)
 {
+  CV_INSTRUMENT_REGION()
+
   kpts.clear();
-  Compute_Determinant_Hessian_Response();
   Find_Scale_Space_Extrema(kpts);
   Do_Subpixel_Refinement(kpts);
 }
 
 /* ************************************************************************* */
-class MultiscaleDerivativesAKAZEInvoker : public ParallelLoopBody
+
+#ifdef HAVE_OPENCL
+static inline bool
+ocl_compute_determinant(const UMat& Lxx, const UMat& Lxy, const UMat& Lyy,
+  UMat& Ldet, float sigma)
+{
+  const int total = Lxx.rows * Lxx.cols;
+  size_t globalSize[] = {(size_t)total};
+
+  ocl::Kernel ker("AKAZE_compute_determinant", ocl::features2d::akaze_oclsrc);
+  if( ker.empty() )
+    return false;
+
+  return ker.args(
+    ocl::KernelArg::PtrReadOnly(Lxx),
+    ocl::KernelArg::PtrReadOnly(Lxy),
+    ocl::KernelArg::PtrReadOnly(Lyy),
+    ocl::KernelArg::PtrWriteOnly(Ldet),
+    sigma, total).run(1, globalSize, 0, true);
+}
+#endif // HAVE_OPENCL
+
+/**
+ * @brief Compute determinant from hessians
+ * @details Compute Ldet by (Lxx.mul(Lyy) - Lxy.mul(Lxy)) * sigma
+ *
+ * @param Lxx spatial derivates
+ * @param Lxy spatial derivates
+ * @param Lyy spatial derivates
+ * @param Ldet output determinant
+ * @param sigma determinant will be scaled by this sigma
+ */
+static inline void compute_determinant(const UMat& Lxx, const UMat& Lxy, const UMat& Lyy,
+  UMat& Ldet, float sigma)
+{
+  CV_INSTRUMENT_REGION()
+
+  Ldet.create(Lxx.size(), Lxx.type());
+
+  CV_OCL_RUN(true, ocl_compute_determinant(Lxx, Lxy, Lyy, Ldet, sigma));
+
+  // output determinant
+  Mat Mxx = Lxx.getMat(ACCESS_READ), Mxy = Lxy.getMat(ACCESS_READ), Myy = Lyy.getMat(ACCESS_READ);
+  Mat Mdet = Ldet.getMat(ACCESS_WRITE);
+  float *lxx = Mxx.ptr<float>();
+  float *lxy = Mxy.ptr<float>();
+  float *lyy = Myy.ptr<float>();
+  float *ldet = Mdet.ptr<float>();
+  const int total = Lxx.cols * Lxx.rows;
+  for (int j = 0; j < total; j++) {
+    ldet[j] = (lxx[j] * lyy[j] - lxy[j] * lxy[j]) * sigma;
+  }
+
+}
+
+class DeterminantHessianResponse : public ParallelLoopBody
 {
 public:
-    explicit MultiscaleDerivativesAKAZEInvoker(std::vector<TEvolution>& ev, const AKAZEOptions& opt)
+    explicit DeterminantHessianResponse(std::vector<Evolution>& ev)
     : evolution_(&ev)
-    , options_(opt)
   {
   }
 
   void operator()(const Range& range) const
   {
-    std::vector<TEvolution>& evolution = *evolution_;
+    UMat Lxx, Lxy, Lyy;
 
     for (int i = range.start; i < range.end; i++)
     {
-      float ratio = (float)fastpow(2, 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_));
+      Evolution &e = (*evolution_)[i];
+
+      // we cannot use cv:Scharr here, because we need to handle also
+      // kernel sizes other than 3, by default we are using 9x9, 5x5 and 7x7
+
+      // compute kernels
+      Mat DxKx, DxKy, DyKx, DyKy;
+      compute_derivative_kernels(DxKx, DxKy, 1, 0, e.sigma_size);
+      compute_derivative_kernels(DyKx, DyKy, 0, 1, e.sigma_size);
+
+      // compute the multiscale derivatives
+      sepFilter2D(e.Lsmooth, e.Lx, CV_32F, DxKx, DxKy);
+      sepFilter2D(e.Lx, Lxx, CV_32F, DxKx, DxKy);
+      sepFilter2D(e.Lx, Lxy, CV_32F, DyKx, DyKy);
+      sepFilter2D(e.Lsmooth, e.Ly, CV_32F, DyKx, DyKy);
+      sepFilter2D(e.Ly, Lyy, CV_32F, DyKx, DyKy);
+
+      // free Lsmooth to same some space in the pyramid, it is not needed anymore
+      e.Lsmooth.release();
+
+      // compute determinant scaled by sigma
+      float sigma_size_quat = (float)(e.sigma_size * e.sigma_size * e.sigma_size * e.sigma_size);
+      compute_determinant(Lxx, Lxy, Lyy, e.Ldet, sigma_size_quat);
     }
   }
 
 private:
-  std::vector<TEvolution>*  evolution_;
-  AKAZEOptions              options_;
+  std::vector<Evolution>*  evolution_;
 };
 
-/* ************************************************************************* */
-/**
- * @brief This method computes the multiscale derivatives for the nonlinear scale space
- */
-void AKAZEFeatures::Compute_Multiscale_Derivatives(void)
-{
-  parallel_for_(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) {
+  CV_INSTRUMENT_REGION()
 
-  // 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);
-      }
-    }
+  if (ocl::useOpenCL()) {
+    DeterminantHessianResponse body (evolution_);
+    body(Range(0, (int)evolution_.size()));
+  } else {
+    parallel_for_(Range(0, (int)evolution_.size()), DeterminantHessianResponse(evolution_));
   }
 }
 
@@ -255,6 +614,7 @@ void AKAZEFeatures::Compute_Determinant_Hessian_Response(void) {
  */
 void AKAZEFeatures::Find_Scale_Space_Extrema(std::vector<KeyPoint>& kpts)
 {
+  CV_INSTRUMENT_REGION()
 
   float value = 0.0;
   float dist = 0.0, ratio = 0.0, smax = 0.0;
@@ -273,16 +633,17 @@ void AKAZEFeatures::Find_Scale_Space_Extrema(std::vector<KeyPoint>& kpts)
   }
 
   for (size_t i = 0; i < evolution_.size(); i++) {
-    float* prev = evolution_[i].Ldet.ptr<float>(0);
-    float* curr = evolution_[i].Ldet.ptr<float>(1);
-    for (int ix = 1; ix < evolution_[i].Ldet.rows - 1; ix++) {
-      float* next = evolution_[i].Ldet.ptr<float>(ix + 1);
+    Mat Ldet = evolution_[i].Mdet;
+    const float* prev = Ldet.ptr<float>(0);
+    const float* curr = Ldet.ptr<float>(1);
+    for (int ix = 1; ix < Ldet.rows - 1; ix++) {
+      const float* next = Ldet.ptr<float>(ix + 1);
 
-      for (int jx = 1; jx < evolution_[i].Ldet.cols - 1; jx++) {
+      for (int jx = 1; jx < Ldet.cols - 1; jx++) {
         is_extremum = false;
         is_repeated = false;
         is_out = false;
-        value = *(evolution_[i].Ldet.ptr<float>(ix)+jx);
+        value = *(Ldet.ptr<float>(ix)+jx);
 
         // Filter the points with the detector threshold
         if (value > options_.dthreshold && value >= options_.min_dthreshold &&
@@ -335,8 +696,8 @@ void AKAZEFeatures::Find_Scale_Space_Extrema(std::vector<KeyPoint>& kpts)
             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) {
+            if (left_x < 0 || right_x >= Ldet.cols ||
+                up_y < 0 || down_y >= Ldet.rows) {
               is_out = true;
             }
 
@@ -394,6 +755,8 @@ void AKAZEFeatures::Find_Scale_Space_Extrema(std::vector<KeyPoint>& kpts)
  */
 void AKAZEFeatures::Do_Subpixel_Refinement(std::vector<KeyPoint>& kpts)
 {
+  CV_INSTRUMENT_REGION()
+
   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;
@@ -405,26 +768,27 @@ void AKAZEFeatures::Do_Subpixel_Refinement(std::vector<KeyPoint>& kpts)
     ratio = (float)fastpow(2, kpts[i].octave);
     x = fRound(kpts[i].pt.x / ratio);
     y = fRound(kpts[i].pt.y / ratio);
+    Mat Ldet = evolution_[kpts[i].class_id].Mdet;
 
     // 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));
+    Dx = (0.5f)*(*(Ldet.ptr<float>(y)+x + 1)
+        - *(Ldet.ptr<float>(y)+x - 1));
+    Dy = (0.5f)*(*(Ldet.ptr<float>(y + 1) + x)
+        - *(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)));
+    Dxx = (*(Ldet.ptr<float>(y)+x + 1)
+        + *(Ldet.ptr<float>(y)+x - 1)
+        - 2.0f*(*(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)));
+    Dyy = (*(Ldet.ptr<float>(y + 1) + x)
+        + *(Ldet.ptr<float>(y - 1) + x)
+        - 2.0f*(*(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)));
+    Dxy = (0.25f)*(*(Ldet.ptr<float>(y + 1) + x + 1)
+        + (*(Ldet.ptr<float>(y - 1) + x - 1)))
+        - (0.25f)*(*(Ldet.ptr<float>(y - 1) + x + 1)
+        + (*(Ldet.ptr<float>(y + 1) + x - 1)));
 
     // Solve the linear system
     A(0, 0) = Dxx;
@@ -459,7 +823,7 @@ void AKAZEFeatures::Do_Subpixel_Refinement(std::vector<KeyPoint>& kpts)
 class SURF_Descriptor_Upright_64_Invoker : public ParallelLoopBody
 {
 public:
-  SURF_Descriptor_Upright_64_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, std::vector<TEvolution>& evolution)
+  SURF_Descriptor_Upright_64_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, std::vector<Evolution>& evolution)
     : keypoints_(&kpts)
     , descriptors_(&desc)
     , evolution_(&evolution)
@@ -479,13 +843,13 @@ public:
 private:
   std::vector<KeyPoint>* keypoints_;
   Mat*                   descriptors_;
-  std::vector<TEvolution>*   evolution_;
+  std::vector<Evolution>*   evolution_;
 };
 
 class SURF_Descriptor_64_Invoker : public ParallelLoopBody
 {
 public:
-  SURF_Descriptor_64_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, std::vector<TEvolution>& evolution)
+  SURF_Descriptor_64_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, std::vector<Evolution>& evolution)
     : keypoints_(&kpts)
     , descriptors_(&desc)
     , evolution_(&evolution)
@@ -496,7 +860,6 @@ public:
   {
     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));
     }
   }
@@ -506,13 +869,13 @@ public:
 private:
   std::vector<KeyPoint>* keypoints_;
   Mat*                   descriptors_;
-  std::vector<TEvolution>*   evolution_;
+  std::vector<Evolution>*   evolution_;
 };
 
 class MSURF_Upright_Descriptor_64_Invoker : public ParallelLoopBody
 {
 public:
-  MSURF_Upright_Descriptor_64_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, std::vector<TEvolution>& evolution)
+  MSURF_Upright_Descriptor_64_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, std::vector<Evolution>& evolution)
     : keypoints_(&kpts)
     , descriptors_(&desc)
     , evolution_(&evolution)
@@ -532,13 +895,13 @@ public:
 private:
   std::vector<KeyPoint>* keypoints_;
   Mat*                   descriptors_;
-  std::vector<TEvolution>*   evolution_;
+  std::vector<Evolution>*   evolution_;
 };
 
 class MSURF_Descriptor_64_Invoker : public ParallelLoopBody
 {
 public:
-  MSURF_Descriptor_64_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, std::vector<TEvolution>& evolution)
+  MSURF_Descriptor_64_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, std::vector<Evolution>& evolution)
     : keypoints_(&kpts)
     , descriptors_(&desc)
     , evolution_(&evolution)
@@ -558,13 +921,13 @@ public:
 private:
   std::vector<KeyPoint>* keypoints_;
   Mat*                   descriptors_;
-  std::vector<TEvolution>*   evolution_;
+  std::vector<Evolution>*   evolution_;
 };
 
 class Upright_MLDB_Full_Descriptor_Invoker : public ParallelLoopBody
 {
 public:
-  Upright_MLDB_Full_Descriptor_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, std::vector<TEvolution>& evolution, AKAZEOptions& options)
+  Upright_MLDB_Full_Descriptor_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, std::vector<Evolution>& evolution, AKAZEOptions& options)
     : keypoints_(&kpts)
     , descriptors_(&desc)
     , evolution_(&evolution)
@@ -585,7 +948,7 @@ public:
 private:
   std::vector<KeyPoint>* keypoints_;
   Mat*                   descriptors_;
-  std::vector<TEvolution>*   evolution_;
+  std::vector<Evolution>*   evolution_;
   AKAZEOptions*              options_;
 };
 
@@ -594,7 +957,7 @@ class Upright_MLDB_Descriptor_Subset_Invoker : public ParallelLoopBody
 public:
   Upright_MLDB_Descriptor_Subset_Invoker(std::vector<KeyPoint>& kpts,
                                          Mat& desc,
-                                         std::vector<TEvolution>& evolution,
+                                         std::vector<Evolution>& evolution,
                                          AKAZEOptions& options,
                                          Mat descriptorSamples,
                                          Mat descriptorBits)
@@ -620,7 +983,7 @@ public:
 private:
   std::vector<KeyPoint>* keypoints_;
   Mat*                   descriptors_;
-  std::vector<TEvolution>*   evolution_;
+  std::vector<Evolution>*   evolution_;
   AKAZEOptions*              options_;
 
   Mat descriptorSamples_;  // List of positions in the grids to sample LDB bits from.
@@ -630,7 +993,7 @@ private:
 class MLDB_Full_Descriptor_Invoker : public ParallelLoopBody
 {
 public:
-  MLDB_Full_Descriptor_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, std::vector<TEvolution>& evolution, AKAZEOptions& options)
+  MLDB_Full_Descriptor_Invoker(std::vector<KeyPoint>& kpts, Mat& desc, std::vector<Evolution>& evolution, AKAZEOptions& options)
     : keypoints_(&kpts)
     , descriptors_(&desc)
     , evolution_(&evolution)
@@ -655,7 +1018,7 @@ public:
 private:
   std::vector<KeyPoint>* keypoints_;
   Mat*                   descriptors_;
-  std::vector<TEvolution>*   evolution_;
+  std::vector<Evolution>*   evolution_;
   AKAZEOptions*              options_;
 };
 
@@ -664,7 +1027,7 @@ class MLDB_Descriptor_Subset_Invoker : public ParallelLoopBody
 public:
   MLDB_Descriptor_Subset_Invoker(std::vector<KeyPoint>& kpts,
                                  Mat& desc,
-                                 std::vector<TEvolution>& evolution,
+                                 std::vector<Evolution>& evolution,
                                  AKAZEOptions& options,
                                  Mat descriptorSamples,
                                  Mat descriptorBits)
@@ -681,7 +1044,6 @@ public:
   {
     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));
     }
   }
@@ -691,7 +1053,7 @@ public:
 private:
   std::vector<KeyPoint>* keypoints_;
   Mat*                   descriptors_;
-  std::vector<TEvolution>*   evolution_;
+  std::vector<Evolution>*   evolution_;
   AKAZEOptions*              options_;
 
   Mat descriptorSamples_;  // List of positions in the grids to sample LDB bits from.
@@ -703,8 +1065,10 @@ private:
  * @param kpts Vector of detected keypoints
  * @param desc Matrix to store the descriptors
  */
-void AKAZEFeatures::Compute_Descriptors(std::vector<KeyPoint>& kpts, Mat& desc)
+void AKAZEFeatures::Compute_Descriptors(std::vector<KeyPoint>& kpts, OutputArray descriptors)
 {
+  CV_INSTRUMENT_REGION()
+
   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()));
@@ -712,20 +1076,22 @@ void AKAZEFeatures::Compute_Descriptors(std::vector<KeyPoint>& kpts, Mat& desc)
 
   // Allocate memory for the matrix with the descriptors
   if (options_.descriptor < AKAZE::DESCRIPTOR_MLDB_UPRIGHT) {
-    desc = Mat::zeros((int)kpts.size(), 64, CV_32FC1);
+    descriptors.create((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 = Mat::zeros((int)kpts.size(), (int)ceil(t / 8.), CV_8UC1);
+      descriptors.create((int)kpts.size(), (int)ceil(t / 8.), CV_8UC1);
     }
     else {
       // We use the random bit selection length binary descriptor
-      desc = Mat::zeros((int)kpts.size(), (int)ceil(options_.descriptor_size / 8.), CV_8UC1);
+      descriptors.create((int)kpts.size(), (int)ceil(options_.descriptor_size / 8.), CV_8UC1);
     }
   }
 
+  Mat desc = descriptors.getMat();
+
   switch (options_.descriptor)
   {
     case AKAZE::DESCRIPTOR_KAZE_UPRIGHT: // Upright descriptors, not invariant to rotation
@@ -759,12 +1125,20 @@ void AKAZEFeatures::Compute_Descriptors(std::vector<KeyPoint>& kpts, Mat& desc)
 
 /* ************************************************************************* */
 /**
- * @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
+ * @brief This function samples the derivative responses Lx and Ly for the points
+ * within the radius of 6*scale from (x0, y0), then multiply 2D Gaussian weight
+ * @param Lx Horizontal derivative
+ * @param Ly Vertical derivative
+ * @param x0 X-coordinate of the center point
+ * @param y0 Y-coordinate of the center point
+ * @param scale The sampling step
+ * @param resX Output array of the weighted horizontal derivative responses
+ * @param resY Output array of the weighted vertical derivative responses
  */
-void AKAZEFeatures::Compute_Main_Orientation(KeyPoint& kpt, const std::vector<TEvolution>& evolution_)
+static inline
+void Sample_Derivative_Response_Radius6(const Mat &Lx, const Mat &Ly,
+                                  const int x0, const int y0, const int scale,
+                                  float *resX, float *resY)
 {
     /* ************************************************************************* */
     /// Lookup table for 2d gaussian (sigma = 2.5) where (0,0) is top left and (6,6) is bottom right
@@ -778,79 +1152,196 @@ void AKAZEFeatures::Compute_Main_Orientation(KeyPoint& kpt, const std::vector<TE
         { 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 }
     };
+    static const int id[] = { 6, 5, 4, 3, 2, 1, 0, 1, 2, 3, 4, 5, 6 };
+    static const struct gtable
+    {
+      float weight[109];
+      int8_t xidx[109];
+      int8_t yidx[109];
 
-  int ix = 0, iy = 0, idx = 0, s = 0, level = 0;
-  float xf = 0.0, yf = 0.0, gweight = 0.0, ratio = 0.0;
-  const int ang_size = 109;
-  float resX[ang_size] = {0}, resY[ang_size] = {0}, Ang[ang_size];
-  const int id[] = { 6, 5, 4, 3, 2, 1, 0, 1, 2, 3, 4, 5, 6 };
+      explicit gtable(void)
+      {
+        // Generate the weight and indices by one-time initialization
+        int k = 0;
+        for (int i = -6; i <= 6; ++i) {
+          for (int j = -6; j <= 6; ++j) {
+            if (i*i + j*j < 36) {
+              weight[k] = gauss25[id[i + 6]][id[j + 6]];
+              yidx[k] = static_cast<int8_t>(i);
+              xidx[k] = static_cast<int8_t>(j);
+              ++k;
+            }
+          }
+        }
+        CV_DbgAssert(k == 109);
+      }
+    } g;
+
+  const float * lx = Lx.ptr<float>(0);
+  const float * ly = Ly.ptr<float>(0);
+  int cols = Lx.cols;
+
+  for (int i = 0; i < 109; i++) {
+    int j = (y0 + g.yidx[i] * scale) * cols + (x0 + g.xidx[i] * scale);
+
+    resX[i] = g.weight[i] * lx[j];
+    resY[i] = g.weight[i] * ly[j];
+
+    CV_DbgAssert(isfinite(resX[i]));
+    CV_DbgAssert(isfinite(resY[i]));
+  }
+}
+
+/**
+ * @brief This function sorts a[] by quantized float values
+ * @param a[] Input floating point array to sort
+ * @param n The length of a[]
+ * @param quantum The interval to convert a[i]'s float values to integers
+ * @param max The upper bound of a[], meaning a[i] must be in [0, max]
+ * @param idx[] Output array of the indices: a[idx[i]] forms a sorted array
+ * @param cum[] Output array of the starting indices of quantized floats
+ * @note The values of a[] in [k*quantum, (k + 1)*quantum) is labeled by
+ * the integer k, which is calculated by floor(a[i]/quantum).  After sorting,
+ * the values from a[idx[cum[k]]] to a[idx[cum[k+1]-1]] are all labeled by k.
+ * This sorting is unstable to reduce the memory access.
+ */
+static inline
+void quantized_counting_sort(const float a[], const int n,
+                             const float quantum, const float max,
+                             uint8_t idx[], uint8_t cum[])
+{
+  const int nkeys = (int)(max / quantum);
+
+  // The size of cum[] must be nkeys + 1
+  memset(cum, 0, nkeys + 1);
 
-  // Variables for computing the dominant direction
-  float sumX = 0.0, sumY = 0.0, max = 0.0, ang1 = 0.0, ang2 = 0.0;
+  // Count up the quantized values
+  for (int i = 0; i < n; i++)
+    cum[(int)(a[i] / quantum)]++;
 
+  // Compute the inclusive prefix sum i.e. the end indices; cum[nkeys] is the total
+  for (int i = 1; i <= nkeys; i++)
+    cum[i] += cum[i - 1];
+
+  // Generate the sorted indices; cum[] becomes the exclusive prefix sum i.e. the start indices of keys
+  for (int i = 0; i < n; i++)
+    idx[--cum[(int)(a[i] / quantum)]] = static_cast<uint8_t>(i);
+}
+
+/**
+ * @brief This function 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
+ */
+static inline
+void Compute_Main_Orientation(KeyPoint& kpt, const std::vector<Evolution>& evolution)
+{
+  // get the right evolution level for this keypoint
+  const Evolution& e = evolution[kpt.class_id];
   // 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;
+  int scale = fRound(0.5f * kpt.size / e.octave_ratio);
+  int x0 = fRound(kpt.pt.x / e.octave_ratio);
+  int y0 = fRound(kpt.pt.y / e.octave_ratio);
+
+  // Sample derivatives responses for the points within radius of 6*scale
+  const int ang_size = 109;
+  float resX[ang_size], resY[ang_size];
+  Sample_Derivative_Response_Radius6(e.Mx, e.My, x0, y0, scale, resX, resY);
+
+  // Compute the angle of each gradient vector
+  float Ang[ang_size];
+  hal::fastAtan2(resY, resX, Ang, ang_size, false);
+
+  // Sort by the angles; angles are labeled by slices of 0.15 radian
+  const int slices = 42;
+  const float ang_step = (float)(2.0 * CV_PI / slices);
+  uint8_t slice[slices + 1];
+  uint8_t sorted_idx[ang_size];
+  quantized_counting_sort(Ang, ang_size, ang_step, (float)(2.0 * CV_PI), sorted_idx, slice);
+
+  // Find the main angle by sliding a window of 7-slice size(=PI/3) around the keypoint
+  const int win = 7;
+
+  float maxX = 0.0f, maxY = 0.0f;
+  for (int i = slice[0]; i < slice[win]; i++) {
+    maxX += resX[sorted_idx[i]];
+    maxY += resY[sorted_idx[i]];
+  }
+  float maxNorm = maxX * maxX + maxY * maxY;
 
-  // 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);
+  for (int sn = 1; sn <= slices - win; sn++) {
 
-        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));
+    if (slice[sn] == slice[sn - 1] && slice[sn + win] == slice[sn + win - 1])
+      continue;  // The contents of the window didn't change; don't repeat the computation
 
-        ++idx;
-      }
+    float sumX = 0.0f, sumY = 0.0f;
+    for (int i = slice[sn]; i < slice[sn + win]; i++) {
+      sumX += resX[sorted_idx[i]];
+      sumY += resY[sorted_idx[i]];
     }
+
+    float norm = sumX * sumX + sumY * sumY;
+    if (norm > maxNorm)
+        maxNorm = norm, maxX = sumX, maxY = sumY;  // Found bigger one; update
   }
-  hal::fastAtan32f(resY, resX, Ang, ang_size, false);
-  // 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 (int 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) * 180.f / static_cast<float>(CV_PI);
+  for (int sn = slices - win + 1; sn < slices; sn++) {
+    int remain = sn + win - slices;
+
+    if (slice[sn] == slice[sn - 1] && slice[remain] == slice[remain - 1])
+      continue;
+
+    float sumX = 0.0f, sumY = 0.0f;
+    for (int i = slice[sn]; i < slice[slices]; i++) {
+      sumX += resX[sorted_idx[i]];
+      sumY += resY[sorted_idx[i]];
     }
+    for (int i = slice[0]; i < slice[remain]; i++) {
+      sumX += resX[sorted_idx[i]];
+      sumY += resY[sorted_idx[i]];
+    }
+
+    float norm = sumX * sumX + sumY * sumY;
+    if (norm > maxNorm)
+        maxNorm = norm, maxX = sumX, maxY = sumY;
   }
+
+  // Store the final result
+  kpt.angle = fastAtan2(maxY, maxX);
 }
 
-/* ************************************************************************* */
+class ComputeKeypointOrientation : public ParallelLoopBody
+{
+public:
+  ComputeKeypointOrientation(std::vector<KeyPoint>& kpts,
+                             const std::vector<Evolution>& evolution)
+    : keypoints_(&kpts)
+    , evolution_(&evolution)
+  {
+  }
+
+  void operator() (const Range& range) const
+  {
+    for (int i = range.start; i < range.end; i++)
+    {
+      Compute_Main_Orientation((*keypoints_)[i], *evolution_);
+    }
+  }
+private:
+  std::vector<KeyPoint>* keypoints_;
+  const std::vector<Evolution>* evolution_;
+};
+
 /**
  * @brief This method computes the main orientation for a given keypoints
  * @param kpts Input keypoints
  */
 void AKAZEFeatures::Compute_Keypoints_Orientation(std::vector<KeyPoint>& kpts) const
 {
-     for(size_t i = 0; i < kpts.size(); i++)
-         Compute_Main_Orientation(kpts[i], evolution_);
+  CV_INSTRUMENT_REGION()
+
+  parallel_for_(Range(0, (int)kpts.size()), ComputeKeypointOrientation(kpts, evolution_));
 }
 
 /* ************************************************************************* */
@@ -871,12 +1362,12 @@ void MSURF_Upright_Descriptor_64_Invoker::Get_MSURF_Upright_Descriptor_64(const
   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;
+  int scale = 0, dsize = 0;
 
   // Subregion centers for the 4x4 gaussian weighting
   float cx = -0.5f, cy = 0.5f;
 
-  const std::vector<TEvolution>& evolution = *evolution_;
+  const std::vector<Evolution>& evolution = *evolution_;
 
   // Set the descriptor size and the sample and pattern sizes
   dsize = 64;
@@ -886,7 +1377,9 @@ void MSURF_Upright_Descriptor_64_Invoker::Get_MSURF_Upright_Descriptor_64(const
   // Get the information from the keypoint
   ratio = (float)(1 << kpt.octave);
   scale = fRound(0.5f*kpt.size / ratio);
-  level = kpt.class_id;
+  const int level = kpt.class_id;
+  Mat Lx = evolution[level].Mx;
+  Mat Ly = evolution[level].My;
   yf = kpt.pt.y / ratio;
   xf = kpt.pt.x / ratio;
 
@@ -929,16 +1422,16 @@ void MSURF_Upright_Descriptor_64_Invoker::Get_MSURF_Upright_Descriptor_64(const
           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 = *(Lx.ptr<float>(y1)+x1);
+          res2 = *(Lx.ptr<float>(y1)+x2);
+          res3 = *(Lx.ptr<float>(y2)+x1);
+          res4 = *(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 = *(Ly.ptr<float>(y1)+x1);
+          res2 = *(Ly.ptr<float>(y1)+x2);
+          res3 = *(Ly.ptr<float>(y2)+x1);
+          res4 = *(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;
@@ -994,12 +1487,12 @@ void MSURF_Descriptor_64_Invoker::Get_MSURF_Descriptor_64(const KeyPoint& kpt, f
   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;
+  int scale = 0, dsize = 0;
 
   // Subregion centers for the 4x4 gaussian weighting
   float cx = -0.5f, cy = 0.5f;
 
-  const std::vector<TEvolution>& evolution = *evolution_;
+  const std::vector<Evolution>& evolution = *evolution_;
 
   // Set the descriptor size and the sample and pattern sizes
   dsize = 64;
@@ -1010,7 +1503,9 @@ void MSURF_Descriptor_64_Invoker::Get_MSURF_Descriptor_64(const KeyPoint& kpt, f
   ratio = (float)(1 << kpt.octave);
   scale = fRound(0.5f*kpt.size / ratio);
   angle = (kpt.angle * static_cast<float>(CV_PI)) / 180.f;
-  level = kpt.class_id;
+  const int level = kpt.class_id;
+  Mat Lx = evolution[level].Mx;
+  Mat Ly = evolution[level].My;
   yf = kpt.pt.y / ratio;
   xf = kpt.pt.x / ratio;
   co = cos(angle);
@@ -1055,26 +1550,26 @@ void MSURF_Descriptor_64_Invoker::Get_MSURF_Descriptor_64(const KeyPoint& kpt, f
 
           // fix crash: indexing with out-of-bounds index, this might happen near the edges of image
           // clip values so they fit into the image
-          const MatSize& size = evolution[level].Lx.size;
+          const MatSize& size = Lx.size;
           y1 = min(max(0, y1), size[0] - 1);
           x1 = min(max(0, x1), size[1] - 1);
           y2 = min(max(0, y2), size[0] - 1);
           x2 = min(max(0, x2), size[1] - 1);
-          CV_DbgAssert(evolution[level].Lx.size == evolution[level].Ly.size);
+          CV_DbgAssert(Lx.size == Ly.size);
 
           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 = *(Lx.ptr<float>(y1, x1));
+          res2 = *(Lx.ptr<float>(y1, x2));
+          res3 = *(Lx.ptr<float>(y2, x1));
+          res4 = *(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 = *(Ly.ptr<float>(y1, x1));
+          res2 = *(Ly.ptr<float>(y1, x2));
+          res3 = *(Ly.ptr<float>(y2, x1));
+          res4 = *(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
@@ -1125,23 +1620,26 @@ void Upright_MLDB_Full_Descriptor_Invoker::Get_Upright_MLDB_Full_Descriptor(cons
   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;
-  int level = 0, nsamples = 0, scale = 0;
+  int nsamples = 0, scale = 0;
   int dcount1 = 0, dcount2 = 0;
 
   const AKAZEOptions & options = *options_;
-  const std::vector<TEvolution>& evolution = *evolution_;
+  const std::vector<Evolution>& evolution = *evolution_;
 
   // Matrices for the M-LDB descriptor
   Mat values[3] = {
-    Mat::zeros(4, options.descriptor_channels, CV_32FC1),
-    Mat::zeros(9, options.descriptor_channels, CV_32FC1),
-    Mat::zeros(16, options.descriptor_channels, CV_32FC1)
+    Mat(4, options.descriptor_channels, CV_32FC1),
+    Mat(9, options.descriptor_channels, CV_32FC1),
+    Mat(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;
+  const int level = kpt.class_id;
+  Mat Lx = evolution[level].Mx;
+  Mat Ly = evolution[level].My;
+  Mat Lt = evolution[level].Mt;
   yf = kpt.pt.y / ratio;
   xf = kpt.pt.x / ratio;
 
@@ -1172,9 +1670,9 @@ void Upright_MLDB_Full_Descriptor_Invoker::Get_Upright_MLDB_Full_Descriptor(cons
             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);
+            ri = *(Lt.ptr<float>(y1)+x1);
+            rx = *(Lx.ptr<float>(y1)+x1);
+            ry = *(Ly.ptr<float>(y1)+x1);
 
             di += ri;
             dx += rx;
@@ -1204,6 +1702,8 @@ void Upright_MLDB_Full_Descriptor_Invoker::Get_Upright_MLDB_Full_Descriptor(cons
         for (int k = 0; k < 3; ++k) {
           if (*(valI + k) > *(valJ + k)) {
             desc[dcount1 / 8] |= (1 << (dcount1 % 8));
+          } else {
+            desc[dcount1 / 8] &= ~(1 << (dcount1 % 8));
           }
           dcount1++;
         }
@@ -1213,13 +1713,16 @@ void Upright_MLDB_Full_Descriptor_Invoker::Get_Upright_MLDB_Full_Descriptor(cons
   } // for (int z = 0; z < 3; z++)
 }
 
-void MLDB_Full_Descriptor_Invoker::MLDB_Fill_Values(float* values, int sample_step, int level,
+void MLDB_Full_Descriptor_Invoker::MLDB_Fill_Values(float* values, int sample_step, const int level,
                                                     float xf, float yf, float co, float si, float scale) const
 {
-    const std::vector<TEvolution>& evolution = *evolution_;
+    const std::vector<Evolution>& evolution = *evolution_;
     int pattern_size = options_->descriptor_pattern_size;
     int chan = options_->descriptor_channels;
     int valpos = 0;
+    Mat Lx = evolution[level].Mx;
+    Mat Ly = evolution[level].My;
+    Mat Lt = evolution[level].Mt;
 
     for (int i = -pattern_size; i < pattern_size; i += sample_step) {
         for (int j = -pattern_size; j < pattern_size; j += sample_step) {
@@ -1237,18 +1740,18 @@ void MLDB_Full_Descriptor_Invoker::MLDB_Fill_Values(float* values, int sample_st
 
                 // fix crash: indexing with out-of-bounds index, this might happen near the edges of image
                 // clip values so they fit into the image
-                const MatSize& size = evolution[level].Lt.size;
-                CV_DbgAssert(size == evolution[level].Lx.size &&
-                             size == evolution[level].Ly.size);
+                const MatSize& size = Lt.size;
+                CV_DbgAssert(size == Lx.size &&
+                             size == Ly.size);
                 y1 = min(max(0, y1), size[0] - 1);
                 x1 = min(max(0, x1), size[1] - 1);
 
-                float ri = *(evolution[level].Lt.ptr<float>(y1, x1));
+                float ri = *(Lt.ptr<float>(y1, x1));
                 di += ri;
 
                 if(chan > 1) {
-                    float rx = *(evolution[level].Lx.ptr<float>(y1, x1));
-                    float ry = *(evolution[level].Ly.ptr<float>(y1, x1));
+                    float rx = *(Lx.ptr<float>(y1, x1));
+                    float ry = *(Ly.ptr<float>(y1, x1));
                     if (chan == 2) {
                       dx += sqrtf(rx*rx + ry*ry);
                     }
@@ -1290,8 +1793,13 @@ void MLDB_Full_Descriptor_Invoker::MLDB_Binary_Comparisons(float* values, unsign
         for (int i = 0; i < count; i++) {
             int ival = ivalues[chan * i + pos];
             for (int j = i + 1; j < count; j++) {
-                int res = ival > ivalues[chan * j + pos];
-                desc[dpos >> 3] |= (res << (dpos & 7));
+                if (ival > ivalues[chan * j + pos]) {
+                  desc[dpos >> 3] |= (1 << (dpos & 7));
+                }
+                else {
+                  desc[dpos >> 3] &= ~(1 << (dpos & 7));
+                }
+
                 dpos++;
             }
         }
@@ -1347,20 +1855,23 @@ void MLDB_Descriptor_Subset_Invoker::Get_MLDB_Descriptor_Subset(const KeyPoint&
   int x1 = 0, y1 = 0;
 
   const AKAZEOptions & options = *options_;
-  const std::vector<TEvolution>& evolution = *evolution_;
+  const std::vector<Evolution>& 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 * static_cast<float>(CV_PI)) / 180.f;
-  int level = kpt.class_id;
+  const int level = kpt.class_id;
+  Mat Lx = evolution[level].Mx;
+  Mat Ly = evolution[level].My;
+  Mat Lt = evolution[level].Mt;
   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
-  Mat values = Mat_<float>::zeros((4 + 9 + 16)*options.descriptor_channels, 1);
+  Mat values((4 + 9 + 16)*options.descriptor_channels, 1, CV_32FC1);
 
   // Sample everything, but only do the comparisons
   vector<int> steps(3);
@@ -1385,11 +1896,11 @@ void MLDB_Descriptor_Subset_Invoker::Get_MLDB_Descriptor_Subset(const KeyPoint&
         y1 = fRound(sample_y);
         x1 = fRound(sample_x);
 
-        di += *(evolution[level].Lt.ptr<float>(y1)+x1);
+        di += *(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);
+          rx = *(Lx.ptr<float>(y1)+x1);
+          ry = *(Ly.ptr<float>(y1)+x1);
 
           if (options.descriptor_channels == 2) {
             dx += sqrtf(rx*rx + ry*ry);
@@ -1421,6 +1932,8 @@ void MLDB_Descriptor_Subset_Invoker::Get_MLDB_Descriptor_Subset(const KeyPoint&
   for (int i = 0; i<descriptorBits_.rows; i++) {
     if (vals[comps[2 * i]] > vals[comps[2 * i + 1]]) {
       desc[i / 8] |= (1 << (i % 8));
+    } else {
+      desc[i / 8] &= ~(1 << (i % 8));
     }
   }
 }
@@ -1441,17 +1954,20 @@ void Upright_MLDB_Descriptor_Subset_Invoker::Get_Upright_MLDB_Descriptor_Subset(
   int x1 = 0, y1 = 0;
 
   const AKAZEOptions & options = *options_;
-  const std::vector<TEvolution>& evolution = *evolution_;
+  const std::vector<Evolution>& 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;
+  const int level = kpt.class_id;
+  Mat Lx = evolution[level].Mx;
+  Mat Ly = evolution[level].My;
+  Mat Lt = evolution[level].Mt;
   float yf = kpt.pt.y / ratio;
   float xf = kpt.pt.x / ratio;
 
   // Allocate memory for the matrix of values
-  Mat values = Mat_<float>::zeros((4 + 9 + 16)*options.descriptor_channels, 1);
+  Mat values ((4 + 9 + 16)*options.descriptor_channels, 1, CV_32FC1);
 
   vector<int> steps(3);
   steps.at(0) = options.descriptor_pattern_size;
@@ -1472,11 +1988,11 @@ void Upright_MLDB_Descriptor_Subset_Invoker::Get_Upright_MLDB_Descriptor_Subset(
 
         y1 = fRound(sample_y);
         x1 = fRound(sample_x);
-        di += *(evolution[level].Lt.ptr<float>(y1)+x1);
+        di += *(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);
+          rx = *(Lx.ptr<float>(y1)+x1);
+          ry = *(Ly.ptr<float>(y1)+x1);
 
           if (options.descriptor_channels == 2) {
             dx += sqrtf(rx*rx + ry*ry);
@@ -1507,6 +2023,8 @@ void Upright_MLDB_Descriptor_Subset_Invoker::Get_Upright_MLDB_Descriptor_Subset(
   for (int i = 0; i<descriptorBits_.rows; i++) {
     if (vals[comps[2 * i]] > vals[comps[2 * i + 1]]) {
       desc[i / 8] |= (1 << (i % 8));
+    } else {
+      desc[i / 8] &= ~(1 << (i % 8));
     }
   }
 }

modules/features2d/src/kaze/AKAZEFeatures.h

@@ -12,11 +12,40 @@
 /* ************************************************************************* */
 // Includes
 #include "AKAZEConfig.h"
-#include "TEvolution.h"
 
 namespace cv
 {
 
+/// A-KAZE nonlinear diffusion filtering evolution
+struct Evolution
+{
+  Evolution() {
+    etime = 0.0f;
+    esigma = 0.0f;
+    octave = 0;
+    sublevel = 0;
+    sigma_size = 0;
+  }
+
+  UMat Lx, Ly;           ///< First order spatial derivatives
+  UMat Lt;               ///< Evolution image
+  UMat Lsmooth;          ///< Smoothed image, used only for computing determinant, released afterwards
+  UMat Ldet;             ///< Detector response
+
+  // the same as above, holding CPU mapping to UMats above
+  Mat Mx, My;
+  Mat Mt;
+  Mat Mdet;
+
+  Size size;                ///< Size of the layer
+  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
+  float octave_ratio;       ///< Scaling ratio of this octave. ratio = 2^octave
+};
+
 /* ************************************************************************* */
 // AKAZE Class Declaration
 class AKAZEFeatures {
@@ -24,7 +53,7 @@ class AKAZEFeatures {
 private:
 
   AKAZEOptions options_;                ///< Configuration options for AKAZE
-  std::vector<TEvolution> evolution_;        ///< Vector of nonlinear diffusion evolution
+  std::vector<Evolution> evolution_;        ///< Vector of nonlinear diffusion evolution
 
   /// FED parameters
   int ncycles_;                  ///< Number of cycles
@@ -44,16 +73,14 @@ public:
 
   /// Scale Space methods
   void Allocate_Memory_Evolution();
-  int Create_Nonlinear_Scale_Space(const cv::Mat& img);
+  void Create_Nonlinear_Scale_Space(InputArray 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_);
+  void Compute_Descriptors(std::vector<cv::KeyPoint>& kpts, OutputArray desc);
   void Compute_Keypoints_Orientation(std::vector<cv::KeyPoint>& kpts) const;
 };
 

modules/features2d/src/kaze/TEvolution.h

@@ -28,6 +28,7 @@ struct TEvolution
   Mat Lt;               ///< Evolution image
   Mat Lsmooth;          ///< Smoothed image
   Mat Ldet;             ///< Detector response
+
   float etime;              ///< Evolution time
   float esigma;             ///< Evolution sigma. For linear diffusion t = sigma^2 / 2
   int octave;               ///< Image octave

modules/features2d/src/kaze/nldiffusion_functions.cpp

@@ -91,7 +91,11 @@ void image_derivatives_scharr(const cv::Mat& src, cv::Mat& dst, int xorder, int
  * @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) {
+void pm_g1(InputArray _Lx, InputArray _Ly, OutputArray _dst, float k) {
+  _dst.create(_Lx.size(), _Lx.type());
+  Mat Lx = _Lx.getMat();
+  Mat Ly = _Ly.getMat();
+  Mat dst = _dst.getMat();
 
   Size sz = Lx.size();
   float inv_k = 1.0f / (k*k);
@@ -118,7 +122,13 @@ void pm_g1(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k) {
  * @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) {
+void pm_g2(InputArray _Lx, InputArray _Ly, OutputArray _dst, float k) {
+    CV_INSTRUMENT_REGION()
+
+    _dst.create(_Lx.size(), _Lx.type());
+    Mat Lx = _Lx.getMat();
+    Mat Ly = _Ly.getMat();
+    Mat dst = _dst.getMat();
 
     Size sz = Lx.size();
     dst.create(sz, Lx.type());
@@ -144,7 +154,11 @@ void pm_g2(const cv::Mat &Lx, const cv::Mat& Ly, cv::Mat& dst, float k) {
  * 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) {
+void weickert_diffusivity(InputArray _Lx, InputArray _Ly, OutputArray _dst, float k) {
+  _dst.create(_Lx.size(), _Lx.type());
+  Mat Lx = _Lx.getMat();
+  Mat Ly = _Ly.getMat();
+  Mat dst = _dst.getMat();
 
   Size sz = Lx.size();
   float inv_k = 1.0f / (k*k);
@@ -177,7 +191,11 @@ void weickert_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, fl
 * 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, float k) {
+void charbonnier_diffusivity(InputArray _Lx, InputArray _Ly, OutputArray _dst, float k) {
+  _dst.create(_Lx.size(), _Lx.type());
+  Mat Lx = _Lx.getMat();
+  Mat Ly = _Ly.getMat();
+  Mat dst = _dst.getMat();
 
   Size sz = Lx.size();
   float inv_k = 1.0f / (k*k);
@@ -209,6 +227,7 @@ void charbonnier_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst,
  * @return k contrast factor
  */
 float compute_k_percentile(const cv::Mat& img, float perc, float gscale, int nbins, int ksize_x, int ksize_y) {
+    CV_INSTRUMENT_REGION()
 
     int nbin = 0, nelements = 0, nthreshold = 0, k = 0;
     float kperc = 0.0, modg = 0.0;
@@ -307,6 +326,7 @@ void compute_scharr_derivatives(const cv::Mat& src, cv::Mat& dst, int xorder, in
  * @param scale_ Scale factor or derivative size
  */
 void compute_derivative_kernels(cv::OutputArray _kx, cv::OutputArray _ky, int dx, int dy, int scale) {
+    CV_INSTRUMENT_REGION()
 
     int ksize = 3 + 2 * (scale - 1);
 
@@ -320,6 +340,7 @@ void compute_derivative_kernels(cv::OutputArray _kx, cv::OutputArray _ky, int dx
     _ky.create(ksize, 1, CV_32F, -1, true);
     Mat kx = _kx.getMat();
     Mat ky = _ky.getMat();
+    std::vector<float> kerI;
 
     float w = 10.0f / 3.0f;
     float norm = 1.0f / (2.0f*scale*(w + 2.0f));
@@ -327,7 +348,7 @@ void compute_derivative_kernels(cv::OutputArray _kx, cv::OutputArray _ky, int dx
     for (int k = 0; k < 2; k++) {
         Mat* kernel = k == 0 ? &kx : &ky;
         int order = k == 0 ? dx : dy;
-        std::vector<float> kerI(ksize, 0.0f);
+        kerI.assign(ksize, 0.0f);
 
         if (order == 0) {
             kerI[0] = norm, kerI[ksize / 2] = w*norm, kerI[ksize - 1] = norm;
@@ -403,6 +424,7 @@ private:
 * 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) {
+    CV_INSTRUMENT_REGION()
 
     cv::parallel_for_(cv::Range(1, Lstep.rows - 1), Nld_Step_Scalar_Invoker(Ld, c, Lstep, stepsize), (double)Ld.total()/(1 << 16));
 
@@ -472,7 +494,6 @@ void nld_step_scalar(cv::Mat& Ld, const cv::Mat& c, cv::Mat& Lstep, float stepsi
 * @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);

modules/features2d/src/kaze/nldiffusion_functions.h

@@ -21,10 +21,10 @@ namespace cv
 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);
-void charbonnier_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k);
+void pm_g1(InputArray Lx, InputArray Ly, OutputArray dst, float k);
+void pm_g2(InputArray Lx, InputArray Ly, OutputArray dst, float k);
+void weickert_diffusivity(InputArray Lx, InputArray Ly, OutputArray dst, float k);
+void charbonnier_diffusivity(InputArray Lx, InputArray Ly, OutputArray dst, float k);
 
 float compute_k_percentile(const cv::Mat& img, float perc, float gscale, int nbins, int ksize_x, int ksize_y);
 

modules/features2d/src/keypoint.cpp

@@ -156,6 +156,8 @@ private:
 
 void KeyPointsFilter::runByPixelsMask( std::vector<KeyPoint>& keypoints, const Mat& mask )
 {
+    CV_INSTRUMENT_REGION()
+
     if( mask.empty() )
         return;
 

modules/features2d/src/opencl/akaze.cl

+// This file is part of OpenCV project.
+// It is subject to the license terms in the LICENSE file found in the top-level directory
+// of this distribution and at http://opencv.org/license.html
+
+
+/**
+ * @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
+ */
+__kernel void
+AKAZE_pm_g2(__global const float* lx, __global const float* ly, __global float* dst,
+    float k, int size)
+{
+    int i = get_global_id(0);
+    // OpenCV plays with dimensions so we need explicit check for this
+    if (!(i < size))
+    {
+        return;
+    }
+
+    const float k2inv = 1.0f / (k * k);
+    dst[i] = 1.0f / (1.0f + ((lx[i] * lx[i] + ly[i] * ly[i]) * k2inv));
+}
+
+__kernel void
+AKAZE_nld_step_scalar(__global const float* lt, int lt_step, int lt_offset, int rows, int cols,
+    __global const float* lf, __global float* dst, float step_size)
+{
+    /* The labeling scheme for this five star stencil:
+        [    a    ]
+        [ -1 c +1 ]
+        [    b    ]
+    */
+    // column-first indexing
+    int i = get_global_id(1);
+    int j = get_global_id(0);
+
+    // OpenCV plays with dimensions so we need explicit check for this
+    if (!(i < rows && j < cols))
+    {
+        return;
+    }
+
+    // get row indexes
+    int a = (i - 1) * cols;
+    int c = (i    ) * cols;
+    int b = (i + 1) * cols;
+    // compute stencil
+    float res = 0.0f;
+    if (i == 0) // first rows
+    {
+        if (j == 0 || j == (cols - 1))
+        {
+            res = 0.0f;
+        } else
+        {
+            res = (lf[c + j] + lf[c + j + 1])*(lt[c + j + 1] - lt[c + j]) +
+                  (lf[c + j] + lf[c + j - 1])*(lt[c + j - 1] - lt[c + j]) +
+                  (lf[c + j] + lf[b + j    ])*(lt[b + j    ] - lt[c + j]);
+        }
+    } else if (i == (rows - 1)) // last row
+    {
+        if (j == 0 || j == (cols - 1))
+        {
+            res = 0.0f;
+        } else
+        {
+            res = (lf[c + j] + lf[c + j + 1])*(lt[c + j + 1] - lt[c + j]) +
+                  (lf[c + j] + lf[c + j - 1])*(lt[c + j - 1] - lt[c + j]) +
+                  (lf[c + j] + lf[a + j    ])*(lt[a + j    ] - lt[c + j]);
+        }
+    } else // inner rows
+    {
+        if (j == 0) // first column
+        {
+            res = (lf[c + 0] + lf[c + 1])*(lt[c + 1] - lt[c + 0]) +
+                  (lf[c + 0] + lf[b + 0])*(lt[b + 0] - lt[c + 0]) +
+                  (lf[c + 0] + lf[a + 0])*(lt[a + 0] - lt[c + 0]);
+        } else if (j == (cols - 1)) // last column
+        {
+            res = (lf[c + j] + lf[c + j - 1])*(lt[c + j - 1] - lt[c + j]) +
+                  (lf[c + j] + lf[b + j    ])*(lt[b + j    ] - lt[c + j]) +
+                  (lf[c + j] + lf[a + j    ])*(lt[a + j    ] - lt[c + j]);
+        } else // inner stencil
+        {
+            res = (lf[c + j] + lf[c + j + 1])*(lt[c + j + 1] - lt[c + j]) +
+                  (lf[c + j] + lf[c + j - 1])*(lt[c + j - 1] - lt[c + j]) +
+                  (lf[c + j] + lf[b + j    ])*(lt[b + j    ] - lt[c + j]) +
+                  (lf[c + j] + lf[a + j    ])*(lt[a + j    ] - lt[c + j]);
+        }
+    }
+
+    dst[c + j] = res * step_size;
+}
+
+/**
+ * @brief Compute determinant from hessians
+ * @details Compute Ldet by (Lxx.mul(Lyy) - Lxy.mul(Lxy)) * sigma
+ *
+ * @param lxx spatial derivates
+ * @param lxy spatial derivates
+ * @param lyy spatial derivates
+ * @param dst output determinant
+ * @param sigma determinant will be scaled by this sigma
+ */
+__kernel void
+AKAZE_compute_determinant(__global const float* lxx, __global const float* lxy, __global const float* lyy,
+    __global float* dst, float sigma, int size)
+{
+    int i = get_global_id(0);
+    // OpenCV plays with dimensions so we need explicit check for this
+    if (!(i < size))
+    {
+        return;
+    }
+
+    dst[i] = (lxx[i] * lyy[i] - lxy[i] * lxy[i]) * sigma;
+}

modules/features2d/test/test_akaze.cpp

+// This file is part of OpenCV project.
+// It is subject to the license terms in the LICENSE file found in the top-level directory
+// of this distribution and at http://opencv.org/license.html
+
+#include "test_precomp.hpp"
+
+using namespace std;
+using namespace cv;
+
+TEST(Features2d_AKAZE, detect_and_compute_split)
+{
+    Mat testImg(100, 100, CV_8U);
+    RNG rng(101);
+    rng.fill(testImg, RNG::UNIFORM, Scalar(0), Scalar(255), true);
+
+    Ptr<Feature2D> ext = AKAZE::create(AKAZE::DESCRIPTOR_MLDB, 0, 3, 0.001f, 1, 1, KAZE::DIFF_PM_G2);
+    vector<KeyPoint> detAndCompKps;
+    Mat desc;
+    ext->detectAndCompute(testImg, noArray(), detAndCompKps, desc);
+
+    vector<KeyPoint> detKps;
+    ext->detect(testImg, detKps);
+
+    ASSERT_EQ(detKps.size(), detAndCompKps.size());
+
+    for(size_t i = 0; i < detKps.size(); i++)
+        ASSERT_EQ(detKps[i].hash(), detAndCompKps[i].hash());
+}
+
+/**
+ * This test is here to guard propagation of NaNs that happens on this image. NaNs are guarded
+ * by debug asserts in AKAZE, which should fire for you if you are lucky.
+ *
+ * This test also reveals problems with uninitialized memory that happens only on this image.
+ * This is very hard to hit and depends a lot on particular allocator. Run this test in valgrind and check
+ * for uninitialized values if you think you are hitting this problem again.
+ */
+TEST(Features2d_AKAZE, uninitialized_and_nans)
+{
+    Mat b1 = imread(cvtest::TS::ptr()->get_data_path() + "../stitching/b1.png");
+    ASSERT_FALSE(b1.empty());
+
+    vector<KeyPoint> keypoints;
+    Mat desc;
+    Ptr<Feature2D> akaze = AKAZE::create();
+    akaze->detectAndCompute(b1, noArray(), keypoints, desc);
+}

modules/features2d/test/test_descriptors_invariance.cpp

@@ -179,7 +179,7 @@ INSTANTIATE_TEST_CASE_P(AKAZE, DescriptorRotationInvariance,
                         Value(IMAGE_TSUKUBA, AKAZE::create(), AKAZE::create(), 0.99f));
 
 INSTANTIATE_TEST_CASE_P(AKAZE_DESCRIPTOR_KAZE, DescriptorRotationInvariance,
-                        Value(IMAGE_TSUKUBA, AKAZE::create(AKAZE::DESCRIPTOR_KAZE), AKAZE::create(AKAZE::DESCRIPTOR_KAZE), 0.002f));
+                        Value(IMAGE_TSUKUBA, AKAZE::create(AKAZE::DESCRIPTOR_KAZE), AKAZE::create(AKAZE::DESCRIPTOR_KAZE), 0.99f));
 
 /*
  * Descriptor's scale invariance check
@@ -189,4 +189,4 @@ INSTANTIATE_TEST_CASE_P(AKAZE, DescriptorScaleInvariance,
                         Value(IMAGE_BIKES, AKAZE::create(), AKAZE::create(), 0.6f));
 
 INSTANTIATE_TEST_CASE_P(AKAZE_DESCRIPTOR_KAZE, DescriptorScaleInvariance,
-                        Value(IMAGE_BIKES, AKAZE::create(AKAZE::DESCRIPTOR_KAZE), AKAZE::create(AKAZE::DESCRIPTOR_KAZE), 0.0004f));
+                        Value(IMAGE_BIKES, AKAZE::create(AKAZE::DESCRIPTOR_KAZE), AKAZE::create(AKAZE::DESCRIPTOR_KAZE), 0.55f));

modules/features2d/test/test_detectors_regression.cpp

@@ -307,24 +307,3 @@ TEST( Features2d_Detector_AKAZE_DESCRIPTOR_KAZE, regression )
     CV_FeatureDetectorTest test( "detector-akaze-with-kaze-desc", AKAZE::create(AKAZE::DESCRIPTOR_KAZE) );
     test.safe_run();
 }
-
-
-TEST( Features2d_Detector_AKAZE, detect_and_compute_split )
-{
-    Mat testImg(100, 100, CV_8U);
-    RNG rng(101);
-    rng.fill(testImg, RNG::UNIFORM, Scalar(0), Scalar(255), true);
-
-    Ptr<Feature2D> ext = AKAZE::create(AKAZE::DESCRIPTOR_MLDB, 0, 3, 0.001f, 1, 1, KAZE::DIFF_PM_G2);
-    vector<KeyPoint> detAndCompKps;
-    Mat desc;
-    ext->detectAndCompute(testImg, noArray(), detAndCompKps, desc);
-
-    vector<KeyPoint> detKps;
-    ext->detect(testImg, detKps);
-
-    ASSERT_EQ(detKps.size(), detAndCompKps.size());
-
-    for(size_t i = 0; i < detKps.size(); i++)
-        ASSERT_EQ(detKps[i].hash(), detAndCompKps[i].hash());
-}

modules/imgproc/src/smooth.cpp

@@ -2386,7 +2386,7 @@ void cv::GaussianBlur( InputArray _src, OutputArray _dst, Size ksize,
     if(sigma1 == 0 && sigma2 == 0 && tegra::useTegra() && tegra::gaussian(src, dst, ksize, borderType))
         return;
 #endif
-    bool useOpenCL = (_dst.isUMat() && _src.dims() <= 2 &&
+    bool useOpenCL = (ocl::useOpenCL() && _dst.isUMat() && _src.dims() <= 2 &&
                ((ksize.width == 3 && ksize.height == 3) ||
                (ksize.width == 5 && ksize.height == 5)) &&
                _src.rows() > ksize.height && _src.cols() > ksize.width);