Added fast_akaze to third party
Signed-off-by: Yash Chainani <yashchainani28@gmail.com>
Change-Id: I7ea2bc5cd3126271f5b04bb8215259044219a675
diff --git a/third_party/akaze/nldiffusion_functions.cpp b/third_party/akaze/nldiffusion_functions.cpp
new file mode 100644
index 0000000..39ec70e
--- /dev/null
+++ b/third_party/akaze/nldiffusion_functions.cpp
@@ -0,0 +1,457 @@
+//=============================================================================
+//
+// 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"
+
+#include <cstdint>
+#include <cstring>
+#include <iostream>
+#include <opencv2/core.hpp>
+#include <opencv2/imgproc.hpp>
+
+// Namespaces
+
+/* ************************************************************************* */
+
+namespace cv {
+using namespace std;
+
+/* ************************************************************************* */
+/**
+ * @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_convolutionV2(const cv::Mat& src, cv::Mat& dst, int ksize_x,
+ int ksize_y, float sigma) {
+ int 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_ = (int)ceil(2.0f * (1.0f + (sigma - 0.8f) / (0.3f)));
+ 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_scharrV2(const cv::Mat& src, cv::Mat& dst, int xorder,
+ int 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_g1V2(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k) {
+ // Compute: dst = exp((Lx.mul(Lx) + Ly.mul(Ly)) / (-k * k))
+
+ const float neg_inv_k2 = -1.0f / (k * k);
+
+ const int total = Lx.rows * Lx.cols;
+ const float* lx = Lx.ptr<float>(0);
+ const float* ly = Ly.ptr<float>(0);
+ float* d = dst.ptr<float>(0);
+
+ for (int i = 0; i < total; i++)
+ d[i] = neg_inv_k2 * (lx[i] * lx[i] + ly[i] * ly[i]);
+
+ exp(dst, 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_g2V2(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k) {
+ // Compute: dst = 1.0f / (1.0f + ((Lx.mul(Lx) + Ly.mul(Ly)) / (k * k)) );
+
+ const float inv_k2 = 1.0f / (k * k);
+
+ const int total = Lx.rows * Lx.cols;
+ const float* lx = Lx.ptr<float>(0);
+ const float* ly = Ly.ptr<float>(0);
+ float* d = dst.ptr<float>(0);
+
+ for (int i = 0; i < total; i++)
+ d[i] = 1.0f / (1.0f + ((lx[i] * lx[i] + ly[i] * ly[i]) * inv_k2));
+}
+
+/* ************************************************************************* */
+/**
+ * @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_diffusivityV2(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst,
+ float k) {
+ // Compute: dst = 1.0f - exp(-3.315f / ((Lx.mul(Lx) + Ly.mul(Ly)) / (k *
+ // k))^4)
+
+ const float inv_k2 = 1.0f / (k * k);
+
+ const int total = Lx.rows * Lx.cols;
+ const float* lx = Lx.ptr<float>(0);
+ const float* ly = Ly.ptr<float>(0);
+ float* d = dst.ptr<float>(0);
+
+ for (int i = 0; i < total; i++) {
+ float dL = inv_k2 * (lx[i] * lx[i] + ly[i] * ly[i]);
+ d[i] = -3.315f / (dL * dL * dL * dL);
+ }
+
+ exp(dst, dst);
+
+ for (int i = 0; i < total; i++) d[i] = 1.0f - d[i];
+}
+
+/* ************************************************************************* */
+/**
+ * @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_diffusivityV2(const cv::Mat& Lx, const cv::Mat& Ly,
+ cv::Mat& dst, float k) {
+ // Compute: dst = 1.0f / sqrt(1.0f + (Lx.mul(Lx) + Ly.mul(Ly)) / (k * k))
+
+ const float inv_k2 = 1.0f / (k * k);
+
+ const int total = Lx.rows * Lx.cols;
+ const float* lx = Lx.ptr<float>(0);
+ const float* ly = Ly.ptr<float>(0);
+ float* d = dst.ptr<float>(0);
+
+ for (int i = 0; i < total; i++)
+ d[i] = 1.0f / sqrtf(1.0f + inv_k2 * (lx[i] * lx[i] + ly[i] * ly[i]));
+}
+
+/* ************************************************************************* */
+/**
+ * @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 perc Percentile of the image gradient histogram (0-1)
+ * @param modgs Temporal vector to hold the gradient norms
+ * @param histogram Temporal vector to hold the gradient histogram
+ * @return k contrast factor
+ */
+float compute_k_percentileV2(const cv::Mat& Lx, const cv::Mat& Ly, float perc,
+ cv::Mat& modgs, cv::Mat& histogram) {
+ const int total = modgs.cols;
+ const int nbins = histogram.cols;
+
+ CV_DbgAssert(total == (Lx.rows - 2) * (Lx.cols - 2));
+ CV_DbgAssert(nbins > 2);
+
+ float* modg = modgs.ptr<float>(0);
+ int32_t* hist = histogram.ptr<int32_t>(0);
+
+ 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++)
+ *modg++ = sqrtf(lx[j] * lx[j] + ly[j] * ly[j]);
+ }
+ modg = modgs.ptr<float>(0);
+
+ // Get the maximum
+ float hmax = 0.0f;
+ for (int i = 0; i < total; i++)
+ if (hmax < modg[i]) hmax = modg[i];
+
+ if (hmax == 0.0f) return 0.03f; // e.g. a blank image
+
+ // Compute the bin numbers: the value range [0, hmax] -> [0, nbins-1]
+ for (int i = 0; i < total; i++) modg[i] *= (nbins - 1) / hmax;
+
+ // Count up
+ std::memset(hist, 0, sizeof(int32_t) * nbins);
+ 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 = nelements + hist[k];
+ }
+
+ return 0.03f;
+}
+
+/* ************************************************************************* */
+/**
+ * @brief Compute Scharr derivative kernels for sizes different than 3
+ * @param _kx Horizontal kernel ues
+ * @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_scharr_derivative_kernelsV2(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, FILTER_SCHARR, 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.0f / 3.0f;
+ float norm = 1.0f / (2.0f * (w + 2.0f));
+
+ std::vector<float> kerI(ksize, 0.0f);
+
+ if (dx == 0) {
+ kerI[0] = norm, kerI[ksize / 2] = w * norm, kerI[ksize - 1] = norm;
+ } else if (dx == 1) {
+ kerI[0] = -1, kerI[ksize / 2] = 0, kerI[ksize - 1] = 1;
+ }
+ Mat(kx.rows, kx.cols, CV_32F, &kerI[0]).copyTo(kx);
+
+ kerI.assign(ksize, 0.0f);
+
+ if (dy == 0) {
+ kerI[0] = norm, kerI[ksize / 2] = w * norm, kerI[ksize - 1] = norm;
+ } else if (dy == 1) {
+ kerI[0] = -1, kerI[ksize / 2] = 0, kerI[ksize - 1] = 1;
+ }
+ Mat(ky.rows, ky.cols, CV_32F, &kerI[0]).copyTo(ky);
+}
+
+inline void nld_step_scalar_one_lane(const cv::Mat& Lt, const cv::Mat& Lf,
+ cv::Mat& Lstep, int idx, int skip) {
+ /* The labeling scheme for this five star stencil:
+ [ a ]
+ [ -1 c +1 ]
+ [ b ]
+ */
+
+ const int cols = Lt.cols - 2;
+ int row = idx;
+
+ const float *lt_a, *lt_c, *lt_b;
+ const float *lf_a, *lf_c, *lf_b;
+ float* dst;
+
+ // 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;
+ dst = Lstep.ptr<float>(0) + 1;
+
+ for (int j = 0; j < cols; j++) {
+ dst[j] = (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]);
+ }
+ row += skip;
+ }
+
+ // Process the middle rows
+ for (; row < Lt.rows - 1; row += skip) {
+ 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
+ dst[0] = (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]);
+
+ lt_a++;
+ lt_c++;
+ lt_b++;
+ lf_a++;
+ lf_c++;
+ lf_b++;
+ dst++;
+
+ // The middle columns
+ for (int j = 0; j < cols; j++) {
+ dst[j] = (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]);
+ }
+
+ // The right-most column
+ dst[cols] = (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]);
+ }
+
+ // Process the bottom row
+ if (row == Lt.rows - 1) {
+ 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;
+ dst = Lstep.ptr<float>(row) + 1;
+
+ for (int j = 0; j < cols; j++) {
+ dst[j] = (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]);
+ }
+ }
+}
+
+/* ************************************************************************* */
+/**
+ * @brief This function computes a scalar non-linear diffusion step
+ * @param Ld Base image in the evolution
+ * @param c Conductivity image
+ * @param Lstep Output image that gives the difference between the current
+ * Ld and the next Ld being evolved
+ * @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_scalarV2(const cv::Mat& Ld, const cv::Mat& c, cv::Mat& Lstep) {
+ nld_step_scalar_one_lane(Ld, c, Lstep, 0, 1);
+}
+
+/* ************************************************************************* */
+/**
+ * @brief This function downsamples the input image using OpenCV resize
+ * @param src Input image to be downsampled
+ * @param dst Output image with half of the resolution of the input image
+ */
+void halfsample_imageV2(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 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_neighbourhoodV2(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;
+}
+
+} // namespace cv
\ No newline at end of file