diff --git a/modules/features2d/src/akaze/AKAZEConfig.h b/modules/features2d/src/akaze/AKAZEConfig.h index 7fed80e2ce..bc3ac93301 100644 --- a/modules/features2d/src/akaze/AKAZEConfig.h +++ b/modules/features2d/src/akaze/AKAZEConfig.h @@ -114,7 +114,7 @@ struct AKAZEOptions { float kcontrast; ///< The contrast factor parameter float kcontrast_percentile; ///< Percentile level for the contrast factor - size_t kcontrast_nbins; ///< Number of bins for the contrast factor histogram + int kcontrast_nbins; ///< Number of bins for the contrast factor histogram bool save_scale_space; ///< Set to true for saving the scale space images bool save_keypoints; ///< Set to true for saving the detected keypoints and descriptors diff --git a/modules/features2d/src/akaze/AKAZEFeatures.cpp b/modules/features2d/src/akaze/AKAZEFeatures.cpp index ece1d7b301..2204f5aba4 100644 --- a/modules/features2d/src/akaze/AKAZEFeatures.cpp +++ b/modules/features2d/src/akaze/AKAZEFeatures.cpp @@ -12,6 +12,7 @@ using namespace std; using namespace cv; +using namespace cv::details::akaze; /* ************************************************************************* */ /** @@ -110,8 +111,7 @@ int AKAZEFeatures::Create_Nonlinear_Scale_Space(const cv::Mat& img) { evolution_[0].Lt.copyTo(evolution_[0].Lsmooth); // First compute the kcontrast factor - options_.kcontrast = compute_k_percentile(img, options_.kcontrast_percentile, - 1.0f, options_.kcontrast_nbins, 0, 0); + options_.kcontrast = compute_k_percentile(img, options_.kcontrast_percentile, 1.0f, options_.kcontrast_nbins, 0, 0); //t2 = cv::getTickCount(); //timing_.kcontrast = 1000.0*(t2 - t1) / cv::getTickFrequency(); diff --git a/modules/features2d/src/akaze/nldiffusion_functions.cpp b/modules/features2d/src/akaze/nldiffusion_functions.cpp index 9ead4ecfe8..e0e2990d29 100644 --- a/modules/features2d/src/akaze/nldiffusion_functions.cpp +++ b/modules/features2d/src/akaze/nldiffusion_functions.cpp @@ -19,368 +19,373 @@ * @author Pablo F. Alcantarilla, Jesus Nuevo */ -#include "nldiffusion_functions.h" +#include "akaze/nldiffusion_functions.h" using namespace std; using namespace cv; -/* ************************************************************************* */ -/** - * @brief This function smoothes an image with a Gaussian kernel - * @param src Input image - * @param dst Output image - * @param ksize_x Kernel size in X-direction (horizontal) - * @param ksize_y Kernel size in Y-direction (vertical) - * @param sigma Kernel standard deviation - */ -void gaussian_2D_convolution(const cv::Mat& src, cv::Mat& dst, const size_t& ksize_x, - const size_t& ksize_y, const float& sigma) { +namespace cv { + namespace details { + namespace akaze { - int ksize_x_ = 0, ksize_y_ = 0; + /* ************************************************************************* */ + /** + * @brief This function smoothes an image with a Gaussian kernel + * @param src Input image + * @param dst Output image + * @param ksize_x Kernel size in X-direction (horizontal) + * @param ksize_y Kernel size in Y-direction (vertical) + * @param sigma Kernel standard deviation + */ + void gaussian_2D_convolution(const cv::Mat& src, cv::Mat& dst, int ksize_x, int ksize_y, float sigma) { - // 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_; - } + int ksize_x_ = 0, ksize_y_ = 0; - // The kernel size must be and odd number - if ((ksize_x_ % 2) == 0) { - ksize_x_ += 1; - } - - if ((ksize_y_ % 2) == 0) { - ksize_y_ += 1; - } - - // Perform the Gaussian Smoothing with border replication - GaussianBlur(src, dst, Size(ksize_x_, ksize_y_), sigma, sigma, BORDER_REPLICATE); -} - -/* ************************************************************************* */ -/** - * @brief This function computes image derivatives with Scharr kernel - * @param src Input image - * @param dst Output image - * @param xorder Derivative order in X-direction (horizontal) - * @param yorder Derivative order in Y-direction (vertical) - * @note Scharr operator approximates better rotation invariance than - * other stencils such as Sobel. See Weickert and Scharr, - * A Scheme for Coherence-Enhancing Diffusion Filtering with Optimized Rotation Invariance, - * Journal of Visual Communication and Image Representation 2002 - */ -void image_derivatives_scharr(const cv::Mat& src, cv::Mat& dst, 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_g1(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k) { - exp(-(Lx.mul(Lx) + Ly.mul(Ly)) / (k*k), dst); -} - -/* ************************************************************************* */ -/** - * @brief This function computes the Perona and Malik conductivity coefficient g2 - * g2 = 1 / (1 + dL^2 / k^2) - * @param Lx First order image derivative in X-direction (horizontal) - * @param Ly First order image derivative in Y-direction (vertical) - * @param dst Output image - * @param k Contrast factor parameter - */ -void pm_g2(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k) { - dst = 1.0 / (1.0 + (Lx.mul(Lx) + Ly.mul(Ly)) / (k*k)); -} - -/* ************************************************************************* */ -/** - * @brief This function computes Weickert conductivity coefficient gw - * @param Lx First order image derivative in X-direction (horizontal) - * @param Ly First order image derivative in Y-direction (vertical) - * @param dst Output image - * @param k Contrast factor parameter - * @note For more information check the following paper: J. Weickert - * Applications of nonlinear diffusion in image processing and computer vision, - * Proceedings of Algorithmy 2000 - */ -void weickert_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k) { - Mat modg; - pow((Lx.mul(Lx) + Ly.mul(Ly)) / (k*k), 4, modg); - cv::exp(-3.315 / modg, dst); - dst = 1.0 - dst; -} - -/* ************************************************************************* */ -/** - * @brief This function computes Charbonnier conductivity coefficient gc - * gc = 1 / sqrt(1 + dL^2 / k^2) - * @param Lx First order image derivative in X-direction (horizontal) - * @param Ly First order image derivative in Y-direction (vertical) - * @param dst Output image - * @param k Contrast factor parameter - * @note For more information check the following paper: J. Weickert - * Applications of nonlinear diffusion in image processing and computer vision, - * Proceedings of Algorithmy 2000 - */ -void charbonnier_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k) { - Mat den; - cv::sqrt(1.0 + (Lx.mul(Lx) + Ly.mul(Ly)) / (k*k), den); - dst = 1.0 / den; -} - -/* ************************************************************************* */ -/** - * @brief This function computes a good empirical value for the k contrast factor - * given an input image, the percentile (0-1), the gradient scale and the number of - * bins in the histogram - * @param img Input image - * @param perc Percentile of the image gradient histogram (0-1) - * @param gscale Scale for computing the image gradient histogram - * @param nbins Number of histogram bins - * @param ksize_x Kernel size in X-direction (horizontal) for the Gaussian smoothing kernel - * @param ksize_y Kernel size in Y-direction (vertical) for the Gaussian smoothing kernel - * @return k contrast factor - */ -float compute_k_percentile(const cv::Mat& img, float perc, float gscale, - size_t nbins, size_t ksize_x, size_t ksize_y) { - - size_t nbin = 0, nelements = 0, nthreshold = 0, k = 0; - float kperc = 0.0, modg = 0.0, lx = 0.0, ly = 0.0; - float npoints = 0.0; - float hmax = 0.0; - - // Create the array for the histogram - std::vector hist(nbins, 0); - - // Create the matrices - cv::Mat gaussian = cv::Mat::zeros(img.rows, img.cols, CV_32F); - cv::Mat Lx = cv::Mat::zeros(img.rows, img.cols, CV_32F); - cv::Mat Ly = cv::Mat::zeros(img.rows, img.cols, CV_32F); - - // Perform the Gaussian convolution - gaussian_2D_convolution(img, gaussian, ksize_x, ksize_y, gscale); - - // Compute the Gaussian derivatives Lx and Ly - image_derivatives_scharr(gaussian, Lx, 1, 0); - image_derivatives_scharr(gaussian, Ly, 0, 1); - - // Skip the borders for computing the histogram - for (int i = 1; i < gaussian.rows - 1; i++) { - for (int j = 1; j < gaussian.cols - 1; j++) { - lx = *(Lx.ptr(i)+j); - ly = *(Ly.ptr(i)+j); - modg = sqrt(lx*lx + ly*ly); - - // Get the maximum - if (modg > hmax) { - hmax = modg; - } - } - } - - // Skip the borders for computing the histogram - for (int i = 1; i < gaussian.rows - 1; i++) { - for (int j = 1; j < gaussian.cols - 1; j++) { - lx = *(Lx.ptr(i)+j); - ly = *(Ly.ptr(i)+j); - modg = sqrt(lx*lx + ly*ly); - - // Find the correspondent bin - if (modg != 0.0) { - nbin = (size_t)floor(nbins*(modg / hmax)); - - if (nbin == nbins) { - nbin--; + // 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_; } - hist[nbin]++; - npoints++; + // 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); } - } - } - // Now find the perc of the histogram percentile - nthreshold = (size_t)(npoints*perc); + /* ************************************************************************* */ + /** + * @brief This function computes image derivatives with Scharr kernel + * @param src Input image + * @param dst Output image + * @param xorder Derivative order in X-direction (horizontal) + * @param yorder Derivative order in Y-direction (vertical) + * @note Scharr operator approximates better rotation invariance than + * other stencils such as Sobel. See Weickert and Scharr, + * A Scheme for Coherence-Enhancing Diffusion Filtering with Optimized Rotation Invariance, + * Journal of Visual Communication and Image Representation 2002 + */ + void image_derivatives_scharr(const cv::Mat& src, cv::Mat& dst, int xorder, int yorder) { + Scharr(src, dst, CV_32F, xorder, yorder, 1.0, 0, BORDER_DEFAULT); + } - for (k = 0; nelements < nthreshold && k < nbins; k++) { - nelements = nelements + hist[k]; - } + /* ************************************************************************* */ + /** + * @brief This function computes the Perona and Malik conductivity coefficient g1 + * g1 = exp(-|dL|^2/k^2) + * @param Lx First order image derivative in X-direction (horizontal) + * @param Ly First order image derivative in Y-direction (vertical) + * @param dst Output image + * @param k Contrast factor parameter + */ + void pm_g1(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k) { + exp(-(Lx.mul(Lx) + Ly.mul(Ly)) / (k*k), dst); + } - if (nelements < nthreshold) { - kperc = 0.03f; - } - else { - kperc = hmax*((float)(k) / (float)nbins); - } + /* ************************************************************************* */ + /** + * @brief This function computes the Perona and Malik conductivity coefficient g2 + * g2 = 1 / (1 + dL^2 / k^2) + * @param Lx First order image derivative in X-direction (horizontal) + * @param Ly First order image derivative in Y-direction (vertical) + * @param dst Output image + * @param k Contrast factor parameter + */ + void pm_g2(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k) { + dst = 1.0 / (1.0 + (Lx.mul(Lx) + Ly.mul(Ly)) / (k*k)); + } - return kperc; -} + /* ************************************************************************* */ + /** + * @brief This function computes Weickert conductivity coefficient gw + * @param Lx First order image derivative in X-direction (horizontal) + * @param Ly First order image derivative in Y-direction (vertical) + * @param dst Output image + * @param k Contrast factor parameter + * @note For more information check the following paper: J. Weickert + * Applications of nonlinear diffusion in image processing and computer vision, + * Proceedings of Algorithmy 2000 + */ + void weickert_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k) { + Mat modg; + pow((Lx.mul(Lx) + Ly.mul(Ly)) / (k*k), 4, modg); + cv::exp(-3.315 / modg, dst); + dst = 1.0 - dst; + } -/* ************************************************************************* */ -/** - * @brief This function computes Scharr image derivatives - * @param src Input image - * @param dst Output image - * @param xorder Derivative order in X-direction (horizontal) - * @param yorder Derivative order in Y-direction (vertical) - * @param scale Scale factor for the derivative size - */ -void compute_scharr_derivatives(const cv::Mat& src, cv::Mat& dst, int xorder, int yorder, int scale) { + /* ************************************************************************* */ + /** + * @brief This function computes Charbonnier conductivity coefficient gc + * gc = 1 / sqrt(1 + dL^2 / k^2) + * @param Lx First order image derivative in X-direction (horizontal) + * @param Ly First order image derivative in Y-direction (vertical) + * @param dst Output image + * @param k Contrast factor parameter + * @note For more information check the following paper: J. Weickert + * Applications of nonlinear diffusion in image processing and computer vision, + * Proceedings of Algorithmy 2000 + */ + void charbonnier_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k) { + Mat den; + cv::sqrt(1.0 + (Lx.mul(Lx) + Ly.mul(Ly)) / (k*k), den); + dst = 1.0 / den; + } - Mat kx, ky; - compute_derivative_kernels(kx, ky, xorder, yorder, scale); - sepFilter2D(src, dst, CV_32F, kx, ky); -} + /* ************************************************************************* */ + /** + * @brief This function computes a good empirical value for the k contrast factor + * given an input image, the percentile (0-1), the gradient scale and the number of + * bins in the histogram + * @param img Input image + * @param perc Percentile of the image gradient histogram (0-1) + * @param gscale Scale for computing the image gradient histogram + * @param nbins Number of histogram bins + * @param ksize_x Kernel size in X-direction (horizontal) for the Gaussian smoothing kernel + * @param ksize_y Kernel size in Y-direction (vertical) for the Gaussian smoothing kernel + * @return k contrast factor + */ + float compute_k_percentile(const cv::Mat& img, float perc, float gscale, int nbins, int ksize_x, int ksize_y) { -/* ************************************************************************* */ -/** - * @brief This function performs a scalar non-linear diffusion step - * @param Ld2 Output image in the evolution - * @param c Conductivity image - * @param Lstep Previous image in the evolution - * @param stepsize The step size in time units - * @note Forward Euler Scheme 3x3 stencil - * The function c is a scalar value that depends on the gradient norm - * dL_by_ds = d(c dL_by_dx)_by_dx + d(c dL_by_dy)_by_dy - */ -void nld_step_scalar(cv::Mat& Ld, const cv::Mat& c, cv::Mat& Lstep, const float& stepsize) { + int nbin = 0, nelements = 0, nthreshold = 0, k = 0; + float kperc = 0.0, modg = 0.0, lx = 0.0, ly = 0.0; + float npoints = 0.0; + float hmax = 0.0; + + // Create the array for the histogram + std::vector hist(nbins, 0); + + // Create the matrices + cv::Mat gaussian = cv::Mat::zeros(img.rows, img.cols, CV_32F); + cv::Mat Lx = cv::Mat::zeros(img.rows, img.cols, CV_32F); + cv::Mat Ly = cv::Mat::zeros(img.rows, img.cols, CV_32F); + + // Perform the Gaussian convolution + gaussian_2D_convolution(img, gaussian, ksize_x, ksize_y, gscale); + + // Compute the Gaussian derivatives Lx and Ly + image_derivatives_scharr(gaussian, Lx, 1, 0); + image_derivatives_scharr(gaussian, Ly, 0, 1); + + // Skip the borders for computing the histogram + for (int i = 1; i < gaussian.rows - 1; i++) { + for (int j = 1; j < gaussian.cols - 1; j++) { + lx = *(Lx.ptr(i)+j); + ly = *(Ly.ptr(i)+j); + modg = sqrt(lx*lx + ly*ly); + + // Get the maximum + if (modg > hmax) { + hmax = modg; + } + } + } + + // Skip the borders for computing the histogram + for (int i = 1; i < gaussian.rows - 1; i++) { + for (int j = 1; j < gaussian.cols - 1; j++) { + lx = *(Lx.ptr(i)+j); + ly = *(Ly.ptr(i)+j); + modg = sqrt(lx*lx + ly*ly); + + // Find the correspondent bin + if (modg != 0.0) { + nbin = (int)floor(nbins*(modg / hmax)); + + if (nbin == nbins) { + nbin--; + } + + hist[nbin]++; + npoints++; + } + } + } + + // Now find the perc of the histogram percentile + nthreshold = (int)(npoints*perc); + + for (k = 0; nelements < nthreshold && k < nbins; k++) { + nelements = nelements + hist[k]; + } + + if (nelements < nthreshold) { + kperc = 0.03f; + } + else { + kperc = hmax*((float)(k) / (float)nbins); + } + + return kperc; + } + + /* ************************************************************************* */ + /** + * @brief This function computes Scharr image derivatives + * @param src Input image + * @param dst Output image + * @param xorder Derivative order in X-direction (horizontal) + * @param yorder Derivative order in Y-direction (vertical) + * @param scale Scale factor for the derivative size + */ + void compute_scharr_derivatives(const cv::Mat& src, cv::Mat& dst, int xorder, int yorder, int scale) { + + Mat kx, ky; + compute_derivative_kernels(kx, ky, xorder, yorder, scale); + sepFilter2D(src, dst, CV_32F, kx, ky); + } + + /* ************************************************************************* */ + /** + * @brief This function performs a scalar non-linear diffusion step + * @param Ld2 Output image in the evolution + * @param c Conductivity image + * @param Lstep Previous image in the evolution + * @param stepsize The step size in time units + * @note Forward Euler Scheme 3x3 stencil + * The function c is a scalar value that depends on the gradient norm + * dL_by_ds = d(c dL_by_dx)_by_dx + d(c dL_by_dy)_by_dy + */ + void nld_step_scalar(cv::Mat& Ld, const cv::Mat& c, cv::Mat& Lstep, const float& stepsize) { #ifdef _OPENMP #pragma omp parallel for schedule(dynamic) #endif - for (int i = 1; i < Lstep.rows - 1; i++) { - for (int j = 1; j < Lstep.cols - 1; j++) { - float xpos = ((*(c.ptr(i)+j)) + (*(c.ptr(i)+j + 1)))*((*(Ld.ptr(i)+j + 1)) - (*(Ld.ptr(i)+j))); - float xneg = ((*(c.ptr(i)+j - 1)) + (*(c.ptr(i)+j)))*((*(Ld.ptr(i)+j)) - (*(Ld.ptr(i)+j - 1))); - float ypos = ((*(c.ptr(i)+j)) + (*(c.ptr(i + 1) + j)))*((*(Ld.ptr(i + 1) + j)) - (*(Ld.ptr(i)+j))); - float yneg = ((*(c.ptr(i - 1) + j)) + (*(c.ptr(i)+j)))*((*(Ld.ptr(i)+j)) - (*(Ld.ptr(i - 1) + j))); - *(Lstep.ptr(i)+j) = 0.5f*stepsize*(xpos - xneg + ypos - yneg); + for (int i = 1; i < Lstep.rows - 1; i++) { + for (int j = 1; j < Lstep.cols - 1; j++) { + float xpos = ((*(c.ptr(i)+j)) + (*(c.ptr(i)+j + 1)))*((*(Ld.ptr(i)+j + 1)) - (*(Ld.ptr(i)+j))); + float xneg = ((*(c.ptr(i)+j - 1)) + (*(c.ptr(i)+j)))*((*(Ld.ptr(i)+j)) - (*(Ld.ptr(i)+j - 1))); + float ypos = ((*(c.ptr(i)+j)) + (*(c.ptr(i + 1) + j)))*((*(Ld.ptr(i + 1) + j)) - (*(Ld.ptr(i)+j))); + float yneg = ((*(c.ptr(i - 1) + j)) + (*(c.ptr(i)+j)))*((*(Ld.ptr(i)+j)) - (*(Ld.ptr(i - 1) + j))); + *(Lstep.ptr(i)+j) = 0.5f*stepsize*(xpos - xneg + ypos - yneg); + } + } + + for (int j = 1; j < Lstep.cols - 1; j++) { + float xpos = ((*(c.ptr(0) + j)) + (*(c.ptr(0) + j + 1)))*((*(Ld.ptr(0) + j + 1)) - (*(Ld.ptr(0) + j))); + float xneg = ((*(c.ptr(0) + j - 1)) + (*(c.ptr(0) + j)))*((*(Ld.ptr(0) + j)) - (*(Ld.ptr(0) + j - 1))); + float ypos = ((*(c.ptr(0) + j)) + (*(c.ptr(1) + j)))*((*(Ld.ptr(1) + j)) - (*(Ld.ptr(0) + j))); + *(Lstep.ptr(0) + j) = 0.5f*stepsize*(xpos - xneg + ypos); + } + + for (int j = 1; j < Lstep.cols - 1; j++) { + float xpos = ((*(c.ptr(Lstep.rows - 1) + j)) + (*(c.ptr(Lstep.rows - 1) + j + 1)))*((*(Ld.ptr(Lstep.rows - 1) + j + 1)) - (*(Ld.ptr(Lstep.rows - 1) + j))); + float xneg = ((*(c.ptr(Lstep.rows - 1) + j - 1)) + (*(c.ptr(Lstep.rows - 1) + j)))*((*(Ld.ptr(Lstep.rows - 1) + j)) - (*(Ld.ptr(Lstep.rows - 1) + j - 1))); + float ypos = ((*(c.ptr(Lstep.rows - 1) + j)) + (*(c.ptr(Lstep.rows - 1) + j)))*((*(Ld.ptr(Lstep.rows - 1) + j)) - (*(Ld.ptr(Lstep.rows - 1) + j))); + float yneg = ((*(c.ptr(Lstep.rows - 2) + j)) + (*(c.ptr(Lstep.rows - 1) + j)))*((*(Ld.ptr(Lstep.rows - 1) + j)) - (*(Ld.ptr(Lstep.rows - 2) + j))); + *(Lstep.ptr(Lstep.rows - 1) + j) = 0.5f*stepsize*(xpos - xneg + ypos - yneg); + } + + for (int i = 1; i < Lstep.rows - 1; i++) { + float xpos = ((*(c.ptr(i))) + (*(c.ptr(i)+1)))*((*(Ld.ptr(i)+1)) - (*(Ld.ptr(i)))); + float xneg = ((*(c.ptr(i))) + (*(c.ptr(i))))*((*(Ld.ptr(i))) - (*(Ld.ptr(i)))); + float ypos = ((*(c.ptr(i))) + (*(c.ptr(i + 1))))*((*(Ld.ptr(i + 1))) - (*(Ld.ptr(i)))); + float yneg = ((*(c.ptr(i - 1))) + (*(c.ptr(i))))*((*(Ld.ptr(i))) - (*(Ld.ptr(i - 1)))); + *(Lstep.ptr(i)) = 0.5f*stepsize*(xpos - xneg + ypos - yneg); + } + + for (int i = 1; i < Lstep.rows - 1; i++) { + float xneg = ((*(c.ptr(i)+Lstep.cols - 2)) + (*(c.ptr(i)+Lstep.cols - 1)))*((*(Ld.ptr(i)+Lstep.cols - 1)) - (*(Ld.ptr(i)+Lstep.cols - 2))); + float ypos = ((*(c.ptr(i)+Lstep.cols - 1)) + (*(c.ptr(i + 1) + Lstep.cols - 1)))*((*(Ld.ptr(i + 1) + Lstep.cols - 1)) - (*(Ld.ptr(i)+Lstep.cols - 1))); + float yneg = ((*(c.ptr(i - 1) + Lstep.cols - 1)) + (*(c.ptr(i)+Lstep.cols - 1)))*((*(Ld.ptr(i)+Lstep.cols - 1)) - (*(Ld.ptr(i - 1) + Lstep.cols - 1))); + *(Lstep.ptr(i)+Lstep.cols - 1) = 0.5f*stepsize*(-xneg + ypos - yneg); + } + + Ld = Ld + Lstep; + } + + /* ************************************************************************* */ + /** + * @brief This function downsamples the input image with the kernel [1/4,1/2,1/4] + * @param img Input image to be downsampled + * @param dst Output image with half of the resolution of the input image + */ + void downsample_image(const cv::Mat& src, cv::Mat& dst) { + + int i1 = 0, j1 = 0, i2 = 0, j2 = 0; + + for (i1 = 1; i1 < src.rows; i1 += 2) { + j2 = 0; + for (j1 = 1; j1 < src.cols; j1 += 2) { + *(dst.ptr(i2)+j2) = 0.5f*(*(src.ptr(i1)+j1)) + 0.25f*(*(src.ptr(i1)+j1 - 1) + *(src.ptr(i1)+j1 + 1)); + j2++; + } + + i2++; + } + } + + /* ************************************************************************* */ + /** + * @brief This function downsamples the input image using OpenCV resize + * @param img Input image to be downsampled + * @param dst Output image with half of the resolution of the input image + */ + void halfsample_image(const cv::Mat& src, cv::Mat& dst) { + + // Make sure the destination image is of the right size + CV_Assert(src.cols / 2 == dst.cols); + CV_Assert(src.rows / 2 == dst.rows); + resize(src, dst, dst.size(), 0, 0, cv::INTER_AREA); + } + + /* ************************************************************************* */ + /** + * @brief Compute Scharr derivative kernels for sizes different than 3 + * @param kx_ The derivative kernel in x-direction + * @param ky_ The derivative kernel in y-direction + * @param dx The derivative order in x-direction + * @param dy The derivative order in y-direction + * @param scale The kernel size + */ + void compute_derivative_kernels(cv::OutputArray kx_, cv::OutputArray ky_, int dx, int dy, int scale) { + + const int ksize = 3 + 2 * (scale - 1); + + // The usual Scharr kernel + if (scale == 1) { + getDerivKernels(kx_, ky_, dx, dy, 0, true, CV_32F); + return; + } + + kx_.create(ksize, 1, CV_32F, -1, true); + ky_.create(ksize, 1, CV_32F, -1, true); + Mat kx = kx_.getMat(); + Mat ky = ky_.getMat(); + + float w = 10.0f / 3.0f; + float norm = 1.0f / (2.0f*scale*(w + 2.0f)); + + for (int k = 0; k < 2; k++) { + Mat* kernel = k == 0 ? &kx : &ky; + int order = k == 0 ? dx : dy; + float kerI[1000]; + + for (int t = 0; t < ksize; t++) { + kerI[t] = 0; + } + + if (order == 0) { + kerI[0] = norm; + kerI[ksize / 2] = w*norm; + kerI[ksize - 1] = norm; + } + else if (order == 1) { + kerI[0] = -1; + kerI[ksize / 2] = 0; + kerI[ksize - 1] = 1; + } + + Mat temp(kernel->rows, kernel->cols, CV_32F, &kerI[0]); + temp.copyTo(*kernel); + } + } } } - - for (int j = 1; j < Lstep.cols - 1; j++) { - float xpos = ((*(c.ptr(0) + j)) + (*(c.ptr(0) + j + 1)))*((*(Ld.ptr(0) + j + 1)) - (*(Ld.ptr(0) + j))); - float xneg = ((*(c.ptr(0) + j - 1)) + (*(c.ptr(0) + j)))*((*(Ld.ptr(0) + j)) - (*(Ld.ptr(0) + j - 1))); - float ypos = ((*(c.ptr(0) + j)) + (*(c.ptr(1) + j)))*((*(Ld.ptr(1) + j)) - (*(Ld.ptr(0) + j))); - *(Lstep.ptr(0) + j) = 0.5f*stepsize*(xpos - xneg + ypos); - } - - for (int j = 1; j < Lstep.cols - 1; j++) { - float xpos = ((*(c.ptr(Lstep.rows - 1) + j)) + (*(c.ptr(Lstep.rows - 1) + j + 1)))*((*(Ld.ptr(Lstep.rows - 1) + j + 1)) - (*(Ld.ptr(Lstep.rows - 1) + j))); - float xneg = ((*(c.ptr(Lstep.rows - 1) + j - 1)) + (*(c.ptr(Lstep.rows - 1) + j)))*((*(Ld.ptr(Lstep.rows - 1) + j)) - (*(Ld.ptr(Lstep.rows - 1) + j - 1))); - float ypos = ((*(c.ptr(Lstep.rows - 1) + j)) + (*(c.ptr(Lstep.rows - 1) + j)))*((*(Ld.ptr(Lstep.rows - 1) + j)) - (*(Ld.ptr(Lstep.rows - 1) + j))); - float yneg = ((*(c.ptr(Lstep.rows - 2) + j)) + (*(c.ptr(Lstep.rows - 1) + j)))*((*(Ld.ptr(Lstep.rows - 1) + j)) - (*(Ld.ptr(Lstep.rows - 2) + j))); - *(Lstep.ptr(Lstep.rows - 1) + j) = 0.5f*stepsize*(xpos - xneg + ypos - yneg); - } - - for (int i = 1; i < Lstep.rows - 1; i++) { - float xpos = ((*(c.ptr(i))) + (*(c.ptr(i)+1)))*((*(Ld.ptr(i)+1)) - (*(Ld.ptr(i)))); - float xneg = ((*(c.ptr(i))) + (*(c.ptr(i))))*((*(Ld.ptr(i))) - (*(Ld.ptr(i)))); - float ypos = ((*(c.ptr(i))) + (*(c.ptr(i + 1))))*((*(Ld.ptr(i + 1))) - (*(Ld.ptr(i)))); - float yneg = ((*(c.ptr(i - 1))) + (*(c.ptr(i))))*((*(Ld.ptr(i))) - (*(Ld.ptr(i - 1)))); - *(Lstep.ptr(i)) = 0.5f*stepsize*(xpos - xneg + ypos - yneg); - } - - for (int i = 1; i < Lstep.rows - 1; i++) { - float xneg = ((*(c.ptr(i)+Lstep.cols - 2)) + (*(c.ptr(i)+Lstep.cols - 1)))*((*(Ld.ptr(i)+Lstep.cols - 1)) - (*(Ld.ptr(i)+Lstep.cols - 2))); - float ypos = ((*(c.ptr(i)+Lstep.cols - 1)) + (*(c.ptr(i + 1) + Lstep.cols - 1)))*((*(Ld.ptr(i + 1) + Lstep.cols - 1)) - (*(Ld.ptr(i)+Lstep.cols - 1))); - float yneg = ((*(c.ptr(i - 1) + Lstep.cols - 1)) + (*(c.ptr(i)+Lstep.cols - 1)))*((*(Ld.ptr(i)+Lstep.cols - 1)) - (*(Ld.ptr(i - 1) + Lstep.cols - 1))); - *(Lstep.ptr(i)+Lstep.cols - 1) = 0.5f*stepsize*(-xneg + ypos - yneg); - } - - Ld = Ld + Lstep; -} - -/* ************************************************************************* */ -/** - * @brief This function downsamples the input image with the kernel [1/4,1/2,1/4] - * @param img Input image to be downsampled - * @param dst Output image with half of the resolution of the input image - */ -void downsample_image(const cv::Mat& src, cv::Mat& dst) { - - int i1 = 0, j1 = 0, i2 = 0, j2 = 0; - - for (i1 = 1; i1 < src.rows; i1 += 2) { - j2 = 0; - for (j1 = 1; j1 < src.cols; j1 += 2) { - *(dst.ptr(i2)+j2) = 0.5f*(*(src.ptr(i1)+j1)) + 0.25f*(*(src.ptr(i1)+j1 - 1) + *(src.ptr(i1)+j1 + 1)); - j2++; - } - - i2++; - } -} - -/* ************************************************************************* */ -/** - * @brief This function downsamples the input image using OpenCV resize - * @param img Input image to be downsampled - * @param dst Output image with half of the resolution of the input image - */ -void halfsample_image(const cv::Mat& src, cv::Mat& dst) { - - // Make sure the destination image is of the right size - CV_Assert(src.cols / 2 == dst.cols); - CV_Assert(src.rows / 2 == dst.rows); - resize(src, dst, dst.size(), 0, 0, cv::INTER_AREA); -} - -/* ************************************************************************* */ -/** - * @brief Compute Scharr derivative kernels for sizes different than 3 - * @param kx_ The derivative kernel in x-direction - * @param ky_ The derivative kernel in y-direction - * @param dx The derivative order in x-direction - * @param dy The derivative order in y-direction - * @param scale The kernel size - */ -void compute_derivative_kernels(cv::OutputArray kx_, cv::OutputArray ky_, int dx, int dy, int scale) { - - const int ksize = 3 + 2 * (scale - 1); - - // The usual Scharr kernel - if (scale == 1) { - getDerivKernels(kx_, ky_, dx, dy, 0, true, CV_32F); - return; - } - - kx_.create(ksize, 1, CV_32F, -1, true); - ky_.create(ksize, 1, CV_32F, -1, true); - Mat kx = kx_.getMat(); - Mat ky = ky_.getMat(); - - float w = 10.0f / 3.0f; - float norm = 1.0f / (2.0f*scale*(w + 2.0f)); - - for (int k = 0; k < 2; k++) { - Mat* kernel = k == 0 ? &kx : &ky; - int order = k == 0 ? dx : dy; - float kerI[1000]; - - for (int t = 0; t < ksize; t++) { - kerI[t] = 0; - } - - if (order == 0) { - kerI[0] = norm; - kerI[ksize / 2] = w*norm; - kerI[ksize - 1] = norm; - } - else if (order == 1) { - kerI[0] = -1; - kerI[ksize / 2] = 0; - kerI[ksize - 1] = 1; - } - - Mat temp(kernel->rows, kernel->cols, CV_32F, &kerI[0]); - temp.copyTo(*kernel); - } -} +} \ No newline at end of file diff --git a/modules/features2d/src/akaze/nldiffusion_functions.h b/modules/features2d/src/akaze/nldiffusion_functions.h index ec0ef2a847..0fab6c59a7 100644 --- a/modules/features2d/src/akaze/nldiffusion_functions.h +++ b/modules/features2d/src/akaze/nldiffusion_functions.h @@ -5,7 +5,8 @@ * @author Pablo F. Alcantarilla, Jesus Nuevo */ -#pragma once +#ifndef AKAZE_NLDIFFUSION_FUNCTIONS_H +#define AKAZE_NLDIFFUSION_FUNCTIONS_H /* ************************************************************************* */ // Includes @@ -13,20 +14,27 @@ /* ************************************************************************* */ // Declaration of functions -void gaussian_2D_convolution(const cv::Mat& src, cv::Mat& dst, const size_t& ksize_x, - const size_t& ksize_y, const float& sigma); -void image_derivatives_scharr(const cv::Mat& src, cv::Mat& dst, - const size_t& xorder, const size_t& yorder); -void pm_g1(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k); -void pm_g2(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k); -void weickert_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k); -void charbonnier_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k); -float compute_k_percentile(const cv::Mat& img, float perc, float gscale, - size_t nbins, size_t ksize_x, size_t ksize_y); -void compute_scharr_derivatives(const cv::Mat& src, cv::Mat& dst, int xorder, int, int scale); -void nld_step_scalar(cv::Mat& Ld, const cv::Mat& c, cv::Mat& Lstep, const float& stepsize); -void downsample_image(const cv::Mat& src, cv::Mat& dst); -void halfsample_image(const cv::Mat& src, cv::Mat& dst); -void compute_derivative_kernels(cv::OutputArray kx_, cv::OutputArray ky_, int dx, int dy, int scale); -bool check_maximum_neighbourhood(const cv::Mat& img, int dsize, float value, - int row, int col, bool same_img); + +namespace cv { + namespace details { + namespace akaze { + + void gaussian_2D_convolution(const cv::Mat& src, cv::Mat& dst, int ksize_x, int ksize_y, float sigma); + void image_derivatives_scharr(const cv::Mat& src, cv::Mat& dst, int xorder, int yorder); + void pm_g1(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k); + void pm_g2(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k); + void weickert_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k); + void charbonnier_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, const float& k); + float compute_k_percentile(const cv::Mat& img, float perc, float gscale, int nbins, int ksize_x, int ksize_y); + void compute_scharr_derivatives(const cv::Mat& src, cv::Mat& dst, int xorder, int, int scale); + void nld_step_scalar(cv::Mat& Ld, const cv::Mat& c, cv::Mat& Lstep, const float& stepsize); + void downsample_image(const cv::Mat& src, cv::Mat& dst); + void halfsample_image(const cv::Mat& src, cv::Mat& dst); + void compute_derivative_kernels(cv::OutputArray kx_, cv::OutputArray ky_, int dx, int dy, int scale); + bool check_maximum_neighbourhood(const cv::Mat& img, int dsize, float value, int row, int col, bool same_img); + + } + } +} + +#endif diff --git a/modules/features2d/src/kaze/KAZEFeatures.cpp b/modules/features2d/src/kaze/KAZEFeatures.cpp index 0fe41aeaa7..4d0127416a 100644 --- a/modules/features2d/src/kaze/KAZEFeatures.cpp +++ b/modules/features2d/src/kaze/KAZEFeatures.cpp @@ -26,6 +26,7 @@ // Namespaces using namespace std; using namespace cv; +using namespace cv::details::kaze; //******************************************************************************* //******************************************************************************* diff --git a/modules/features2d/src/kaze/nldiffusion_functions.cpp b/modules/features2d/src/kaze/nldiffusion_functions.cpp index c2c46d2b7a..23ffaf1f34 100644 --- a/modules/features2d/src/kaze/nldiffusion_functions.cpp +++ b/modules/features2d/src/kaze/nldiffusion_functions.cpp @@ -28,349 +28,355 @@ // Namespaces using namespace std; using namespace cv; +using namespace cv::details::kaze; //************************************************************************************* //************************************************************************************* -/** - * @brief This function smoothes an image with a Gaussian kernel - * @param src Input image - * @param dst Output image - * @param ksize_x Kernel size in X-direction (horizontal) - * @param ksize_y Kernel size in Y-direction (vertical) - * @param sigma Kernel standard deviation - */ -void gaussian_2D_convolution(const cv::Mat& src, cv::Mat& dst, - int ksize_x, int ksize_y, float sigma) { +namespace cv { + namespace details { + namespace kaze { + /** + * @brief This function smoothes an image with a Gaussian kernel + * @param src Input image + * @param dst Output image + * @param ksize_x Kernel size in X-direction (horizontal) + * @param ksize_y Kernel size in Y-direction (vertical) + * @param sigma Kernel standard deviation + */ + void gaussian_2D_convolution(const cv::Mat& src, cv::Mat& dst, + int ksize_x, int ksize_y, float sigma) { - int ksize_x_ = 0, ksize_y_ = 0; + 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_; - } + // 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; - } + // The kernel size must be and odd number + if ((ksize_x_ % 2) == 0) { + ksize_x_ += 1; + } - if ((ksize_y_ % 2) == 0) { - ksize_y_ += 1; - } + if ((ksize_y_ % 2) == 0) { + ksize_y_ += 1; + } - // Perform the Gaussian Smoothing with border replication - GaussianBlur(src,dst,Size(ksize_x_,ksize_y_),sigma,sigma,cv::BORDER_REPLICATE); -} + // Perform the Gaussian Smoothing with border replication + GaussianBlur(src, dst, Size(ksize_x_, ksize_y_), sigma, sigma, cv::BORDER_REPLICATE); + } -//************************************************************************************* -//************************************************************************************* + //************************************************************************************* + //************************************************************************************* -/** - * @brief This function computes the Perona and Malik conductivity coefficient g1 - * g1 = exp(-|dL|^2/k^2) - * @param Lx First order image derivative in X-direction (horizontal) - * @param Ly First order image derivative in Y-direction (vertical) - * @param dst Output image - * @param k Contrast factor parameter - */ -void pm_g1(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k) { - cv::exp(-(Lx.mul(Lx) + Ly.mul(Ly))/(k*k),dst); -} + /** + * @brief This function computes the Perona and Malik conductivity coefficient g1 + * g1 = exp(-|dL|^2/k^2) + * @param Lx First order image derivative in X-direction (horizontal) + * @param Ly First order image derivative in Y-direction (vertical) + * @param dst Output image + * @param k Contrast factor parameter + */ + void pm_g1(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k) { + cv::exp(-(Lx.mul(Lx) + Ly.mul(Ly)) / (k*k), dst); + } -//************************************************************************************* -//************************************************************************************* + //************************************************************************************* + //************************************************************************************* -/** - * @brief This function computes the Perona and Malik conductivity coefficient g2 - * g2 = 1 / (1 + dL^2 / k^2) - * @param Lx First order image derivative in X-direction (horizontal) - * @param Ly First order image derivative in Y-direction (vertical) - * @param dst Output image - * @param k Contrast factor parameter - */ -void pm_g2(const cv::Mat &Lx, const cv::Mat& Ly, cv::Mat& dst, float k) { - dst = 1./(1. + (Lx.mul(Lx) + Ly.mul(Ly))/(k*k)); -} + /** + * @brief This function computes the Perona and Malik conductivity coefficient g2 + * g2 = 1 / (1 + dL^2 / k^2) + * @param Lx First order image derivative in X-direction (horizontal) + * @param Ly First order image derivative in Y-direction (vertical) + * @param dst Output image + * @param k Contrast factor parameter + */ + void pm_g2(const cv::Mat &Lx, const cv::Mat& Ly, cv::Mat& dst, float k) { + dst = 1. / (1. + (Lx.mul(Lx) + Ly.mul(Ly)) / (k*k)); + } -//************************************************************************************* -//************************************************************************************* + //************************************************************************************* + //************************************************************************************* -/** - * @brief This function computes Weickert conductivity coefficient g3 - * @param Lx First order image derivative in X-direction (horizontal) - * @param Ly First order image derivative in Y-direction (vertical) - * @param dst Output image - * @param k Contrast factor parameter - * @note For more information check the following paper: J. Weickert - * Applications of nonlinear diffusion in image processing and computer vision, - * Proceedings of Algorithmy 2000 - */ -void weickert_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k) { - Mat modg; - cv::pow((Lx.mul(Lx) + Ly.mul(Ly))/(k*k),4,modg); - cv::exp(-3.315/modg, dst); - dst = 1.0f - dst; -} + /** + * @brief This function computes Weickert conductivity coefficient g3 + * @param Lx First order image derivative in X-direction (horizontal) + * @param Ly First order image derivative in Y-direction (vertical) + * @param dst Output image + * @param k Contrast factor parameter + * @note For more information check the following paper: J. Weickert + * Applications of nonlinear diffusion in image processing and computer vision, + * Proceedings of Algorithmy 2000 + */ + void weickert_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k) { + Mat modg; + cv::pow((Lx.mul(Lx) + Ly.mul(Ly)) / (k*k), 4, modg); + cv::exp(-3.315 / modg, dst); + dst = 1.0f - dst; + } -//************************************************************************************* -//************************************************************************************* + //************************************************************************************* + //************************************************************************************* -/** - * @brief This function computes a good empirical value for the k contrast factor - * given an input image, the percentile (0-1), the gradient scale and the number of - * bins in the histogram - * @param img Input image - * @param perc Percentile of the image gradient histogram (0-1) - * @param gscale Scale for computing the image gradient histogram - * @param nbins Number of histogram bins - * @param ksize_x Kernel size in X-direction (horizontal) for the Gaussian smoothing kernel - * @param ksize_y Kernel size in Y-direction (vertical) for the Gaussian smoothing kernel - * @return k contrast factor - */ -float compute_k_percentile(const cv::Mat& img, float perc, float gscale, - int nbins, int ksize_x, int ksize_y) { + /** + * @brief This function computes a good empirical value for the k contrast factor + * given an input image, the percentile (0-1), the gradient scale and the number of + * bins in the histogram + * @param img Input image + * @param perc Percentile of the image gradient histogram (0-1) + * @param gscale Scale for computing the image gradient histogram + * @param nbins Number of histogram bins + * @param ksize_x Kernel size in X-direction (horizontal) for the Gaussian smoothing kernel + * @param ksize_y Kernel size in Y-direction (vertical) for the Gaussian smoothing kernel + * @return k contrast factor + */ + float compute_k_percentile(const cv::Mat& img, float perc, float gscale, int nbins, int ksize_x, int ksize_y) { - int nbin = 0, nelements = 0, nthreshold = 0, k = 0; - float kperc = 0.0, modg = 0.0, lx = 0.0, ly = 0.0; - float npoints = 0.0; - float hmax = 0.0; + int nbin = 0, nelements = 0, nthreshold = 0, k = 0; + float kperc = 0.0, modg = 0.0, lx = 0.0, ly = 0.0; + float npoints = 0.0; + float hmax = 0.0; - // Create the array for the histogram - std::vector hist(nbins, 0); + // Create the array for the histogram + std::vector hist(nbins, 0); - // Create the matrices - Mat gaussian = Mat::zeros(img.rows,img.cols,CV_32F); - Mat Lx = Mat::zeros(img.rows,img.cols,CV_32F); - Mat Ly = Mat::zeros(img.rows,img.cols,CV_32F); + // Create the matrices + Mat gaussian = Mat::zeros(img.rows, img.cols, CV_32F); + Mat Lx = Mat::zeros(img.rows, img.cols, CV_32F); + Mat Ly = Mat::zeros(img.rows, img.cols, CV_32F); - // Perform the Gaussian convolution - gaussian_2D_convolution(img,gaussian,ksize_x,ksize_y,gscale); + // Perform the Gaussian convolution + gaussian_2D_convolution(img, gaussian, ksize_x, ksize_y, gscale); - // Compute the Gaussian derivatives Lx and Ly - Scharr(gaussian,Lx,CV_32F,1,0,1,0,cv::BORDER_DEFAULT); - Scharr(gaussian,Ly,CV_32F,0,1,1,0,cv::BORDER_DEFAULT); + // Compute the Gaussian derivatives Lx and Ly + Scharr(gaussian, Lx, CV_32F, 1, 0, 1, 0, cv::BORDER_DEFAULT); + Scharr(gaussian, Ly, CV_32F, 0, 1, 1, 0, cv::BORDER_DEFAULT); - // Skip the borders for computing the histogram - for (int i = 1; i < gaussian.rows-1; i++) { - for (int j = 1; j < gaussian.cols-1; j++) { - lx = *(Lx.ptr(i)+j); - ly = *(Ly.ptr(i)+j); - modg = sqrt(lx*lx + ly*ly); + // Skip the borders for computing the histogram + for (int i = 1; i < gaussian.rows - 1; i++) { + for (int j = 1; j < gaussian.cols - 1; j++) { + lx = *(Lx.ptr(i)+j); + ly = *(Ly.ptr(i)+j); + modg = sqrt(lx*lx + ly*ly); - // Get the maximum - if (modg > hmax) { - hmax = modg; - } - } - } + // Get the maximum + if (modg > hmax) { + hmax = modg; + } + } + } - // Skip the borders for computing the histogram - for (int i = 1; i < gaussian.rows-1; i++) { - for (int j = 1; j < gaussian.cols-1; j++) { - lx = *(Lx.ptr(i)+j); - ly = *(Ly.ptr(i)+j); - modg = sqrt(lx*lx + ly*ly); + // Skip the borders for computing the histogram + for (int i = 1; i < gaussian.rows - 1; i++) { + for (int j = 1; j < gaussian.cols - 1; j++) { + lx = *(Lx.ptr(i)+j); + ly = *(Ly.ptr(i)+j); + modg = sqrt(lx*lx + ly*ly); - // Find the correspondent bin - if (modg != 0.0) { - nbin = (int)floor(nbins*(modg/hmax)); + // Find the correspondent bin + if (modg != 0.0) { + nbin = (int)floor(nbins*(modg / hmax)); - if (nbin == nbins) { - nbin--; - } + if (nbin == nbins) { + nbin--; + } - hist[nbin]++; - npoints++; - } - } - } + hist[nbin]++; + npoints++; + } + } + } - // Now find the perc of the histogram percentile - nthreshold = (size_t)(npoints*perc); + // Now find the perc of the histogram percentile + nthreshold = (size_t)(npoints*perc); - for (k = 0; nelements < nthreshold && k < nbins; k++) { - nelements = nelements + hist[k]; - } + for (k = 0; nelements < nthreshold && k < nbins; k++) { + nelements = nelements + hist[k]; + } - if (nelements < nthreshold) { - kperc = 0.03f; - } - else { - kperc = hmax*((float)(k)/(float)nbins); - } + if (nelements < nthreshold) { + kperc = 0.03f; + } + else { + kperc = hmax*((float)(k) / (float)nbins); + } - return kperc; -} + return kperc; + } -//************************************************************************************* -//************************************************************************************* + //************************************************************************************* + //************************************************************************************* -/** - * @brief This function computes Scharr image derivatives - * @param src Input image - * @param dst Output image - * @param xorder Derivative order in X-direction (horizontal) - * @param yorder Derivative order in Y-direction (vertical) - * @param scale Scale factor or derivative size - */ -void compute_scharr_derivatives(const cv::Mat& src, cv::Mat& dst, - int xorder, int yorder, int scale) { - Mat kx, ky; - compute_derivative_kernels(kx,ky,xorder,yorder,scale); - sepFilter2D(src,dst,CV_32F,kx,ky); -} + /** + * @brief This function computes Scharr image derivatives + * @param src Input image + * @param dst Output image + * @param xorder Derivative order in X-direction (horizontal) + * @param yorder Derivative order in Y-direction (vertical) + * @param scale Scale factor or derivative size + */ + void compute_scharr_derivatives(const cv::Mat& src, cv::Mat& dst, + int xorder, int yorder, int scale) { + Mat kx, ky; + compute_derivative_kernels(kx, ky, xorder, yorder, scale); + sepFilter2D(src, dst, CV_32F, kx, ky); + } -//************************************************************************************* -//************************************************************************************* + //************************************************************************************* + //************************************************************************************* -/** - * @brief Compute derivative kernels for sizes different than 3 - * @param _kx Horizontal kernel values - * @param _ky Vertical kernel values - * @param dx Derivative order in X-direction (horizontal) - * @param dy Derivative order in Y-direction (vertical) - * @param scale_ Scale factor or derivative size - */ -void compute_derivative_kernels(cv::OutputArray _kx, cv::OutputArray _ky, - int dx, int dy, int scale) { + /** + * @brief Compute derivative kernels for sizes different than 3 + * @param _kx Horizontal kernel values + * @param _ky Vertical kernel values + * @param dx Derivative order in X-direction (horizontal) + * @param dy Derivative order in Y-direction (vertical) + * @param scale_ Scale factor or derivative size + */ + void compute_derivative_kernels(cv::OutputArray _kx, cv::OutputArray _ky, + int dx, int dy, int scale) { - int ksize = 3 + 2*(scale-1); + int ksize = 3 + 2 * (scale - 1); - // The standard Scharr kernel - if (scale == 1) { - getDerivKernels(_kx,_ky,dx,dy,0,true,CV_32F); - return; - } + // The standard Scharr kernel + if (scale == 1) { + getDerivKernels(_kx, _ky, dx, dy, 0, true, CV_32F); + return; + } - _kx.create(ksize,1,CV_32F,-1,true); - _ky.create(ksize,1,CV_32F,-1,true); - Mat kx = _kx.getMat(); - Mat ky = _ky.getMat(); + _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*scale*(w+2.0f)); + float w = 10.0f / 3.0f; + float norm = 1.0f / (2.0f*scale*(w + 2.0f)); - for (int k = 0; k < 2; k++) { - Mat* kernel = k == 0 ? &kx : &ky; - int order = k == 0 ? dx : dy; - std::vector kerI(ksize, 0.0f); + for (int k = 0; k < 2; k++) { + Mat* kernel = k == 0 ? &kx : &ky; + int order = k == 0 ? dx : dy; + std::vector kerI(ksize, 0.0f); - if (order == 0) { - kerI[0] = norm, kerI[ksize/2] = w*norm, kerI[ksize-1] = norm; - } - else if (order == 1) { - kerI[0] = -1, kerI[ksize/2] = 0, kerI[ksize-1] = 1; - } + if (order == 0) { + kerI[0] = norm, kerI[ksize / 2] = w*norm, kerI[ksize - 1] = norm; + } + else if (order == 1) { + kerI[0] = -1, kerI[ksize / 2] = 0, kerI[ksize - 1] = 1; + } - Mat temp(kernel->rows,kernel->cols,CV_32F,&kerI[0]); - temp.copyTo(*kernel); - } -} + Mat temp(kernel->rows, kernel->cols, CV_32F, &kerI[0]); + temp.copyTo(*kernel); + } + } -//************************************************************************************* -//************************************************************************************* + //************************************************************************************* + //************************************************************************************* -/** - * @brief This function performs a scalar non-linear diffusion step - * @param Ld2 Output image in the evolution - * @param c Conductivity image - * @param Lstep Previous image in the evolution - * @param stepsize The step size in time units - * @note Forward Euler Scheme 3x3 stencil - * The function c is a scalar value that depends on the gradient norm - * dL_by_ds = d(c dL_by_dx)_by_dx + d(c dL_by_dy)_by_dy - */ -void nld_step_scalar(cv::Mat& Ld, const cv::Mat& c, cv::Mat& Lstep, float stepsize) { + /** + * @brief This function performs a scalar non-linear diffusion step + * @param Ld2 Output image in the evolution + * @param c Conductivity image + * @param Lstep Previous image in the evolution + * @param stepsize The step size in time units + * @note Forward Euler Scheme 3x3 stencil + * The function c is a scalar value that depends on the gradient norm + * dL_by_ds = d(c dL_by_dx)_by_dx + d(c dL_by_dy)_by_dy + */ + void nld_step_scalar(cv::Mat& Ld, const cv::Mat& c, cv::Mat& Lstep, float stepsize) { #ifdef _OPENMP #pragma omp parallel for schedule(dynamic) #endif - for (int i = 1; i < Lstep.rows-1; i++) { - for (int j = 1; j < Lstep.cols-1; j++) { - float xpos = ((*(c.ptr(i)+j))+(*(c.ptr(i)+j+1)))*((*(Ld.ptr(i)+j+1))-(*(Ld.ptr(i)+j))); - float xneg = ((*(c.ptr(i)+j-1))+(*(c.ptr(i)+j)))*((*(Ld.ptr(i)+j))-(*(Ld.ptr(i)+j-1))); - float ypos = ((*(c.ptr(i)+j))+(*(c.ptr(i+1)+j)))*((*(Ld.ptr(i+1)+j))-(*(Ld.ptr(i)+j))); - float yneg = ((*(c.ptr(i-1)+j))+(*(c.ptr(i)+j)))*((*(Ld.ptr(i)+j))-(*(Ld.ptr(i-1)+j))); - *(Lstep.ptr(i)+j) = 0.5f*stepsize*(xpos-xneg + ypos-yneg); - } - } + for (int i = 1; i < Lstep.rows - 1; i++) { + for (int j = 1; j < Lstep.cols - 1; j++) { + float xpos = ((*(c.ptr(i)+j)) + (*(c.ptr(i)+j + 1)))*((*(Ld.ptr(i)+j + 1)) - (*(Ld.ptr(i)+j))); + float xneg = ((*(c.ptr(i)+j - 1)) + (*(c.ptr(i)+j)))*((*(Ld.ptr(i)+j)) - (*(Ld.ptr(i)+j - 1))); + float ypos = ((*(c.ptr(i)+j)) + (*(c.ptr(i + 1) + j)))*((*(Ld.ptr(i + 1) + j)) - (*(Ld.ptr(i)+j))); + float yneg = ((*(c.ptr(i - 1) + j)) + (*(c.ptr(i)+j)))*((*(Ld.ptr(i)+j)) - (*(Ld.ptr(i - 1) + j))); + *(Lstep.ptr(i)+j) = 0.5f*stepsize*(xpos - xneg + ypos - yneg); + } + } - for (int j = 1; j < Lstep.cols-1; j++) { - float xpos = ((*(c.ptr(0)+j))+(*(c.ptr(0)+j+1)))*((*(Ld.ptr(0)+j+1))-(*(Ld.ptr(0)+j))); - float xneg = ((*(c.ptr(0)+j-1))+(*(c.ptr(0)+j)))*((*(Ld.ptr(0)+j))-(*(Ld.ptr(0)+j-1))); - float ypos = ((*(c.ptr(0)+j))+(*(c.ptr(1)+j)))*((*(Ld.ptr(1)+j))-(*(Ld.ptr(0)+j))); - float yneg = ((*(c.ptr(0)+j))+(*(c.ptr(0)+j)))*((*(Ld.ptr(0)+j))-(*(Ld.ptr(0)+j))); - *(Lstep.ptr(0)+j) = 0.5f*stepsize*(xpos-xneg + ypos-yneg); - } + for (int j = 1; j < Lstep.cols - 1; j++) { + float xpos = ((*(c.ptr(0) + j)) + (*(c.ptr(0) + j + 1)))*((*(Ld.ptr(0) + j + 1)) - (*(Ld.ptr(0) + j))); + float xneg = ((*(c.ptr(0) + j - 1)) + (*(c.ptr(0) + j)))*((*(Ld.ptr(0) + j)) - (*(Ld.ptr(0) + j - 1))); + float ypos = ((*(c.ptr(0) + j)) + (*(c.ptr(1) + j)))*((*(Ld.ptr(1) + j)) - (*(Ld.ptr(0) + j))); + float yneg = ((*(c.ptr(0) + j)) + (*(c.ptr(0) + j)))*((*(Ld.ptr(0) + j)) - (*(Ld.ptr(0) + j))); + *(Lstep.ptr(0) + j) = 0.5f*stepsize*(xpos - xneg + ypos - yneg); + } - for (int j = 1; j < Lstep.cols-1; j++) { - float xpos = ((*(c.ptr(Lstep.rows-1)+j))+(*(c.ptr(Lstep.rows-1)+j+1)))*((*(Ld.ptr(Lstep.rows-1)+j+1))-(*(Ld.ptr(Lstep.rows-1)+j))); - float xneg = ((*(c.ptr(Lstep.rows-1)+j-1))+(*(c.ptr(Lstep.rows-1)+j)))*((*(Ld.ptr(Lstep.rows-1)+j))-(*(Ld.ptr(Lstep.rows-1)+j-1))); - float ypos = ((*(c.ptr(Lstep.rows-1)+j))+(*(c.ptr(Lstep.rows-1)+j)))*((*(Ld.ptr(Lstep.rows-1)+j))-(*(Ld.ptr(Lstep.rows-1)+j))); - float yneg = ((*(c.ptr(Lstep.rows-2)+j))+(*(c.ptr(Lstep.rows-1)+j)))*((*(Ld.ptr(Lstep.rows-1)+j))-(*(Ld.ptr(Lstep.rows-2)+j))); - *(Lstep.ptr(Lstep.rows-1)+j) = 0.5f*stepsize*(xpos-xneg + ypos-yneg); - } + for (int j = 1; j < Lstep.cols - 1; j++) { + float xpos = ((*(c.ptr(Lstep.rows - 1) + j)) + (*(c.ptr(Lstep.rows - 1) + j + 1)))*((*(Ld.ptr(Lstep.rows - 1) + j + 1)) - (*(Ld.ptr(Lstep.rows - 1) + j))); + float xneg = ((*(c.ptr(Lstep.rows - 1) + j - 1)) + (*(c.ptr(Lstep.rows - 1) + j)))*((*(Ld.ptr(Lstep.rows - 1) + j)) - (*(Ld.ptr(Lstep.rows - 1) + j - 1))); + float ypos = ((*(c.ptr(Lstep.rows - 1) + j)) + (*(c.ptr(Lstep.rows - 1) + j)))*((*(Ld.ptr(Lstep.rows - 1) + j)) - (*(Ld.ptr(Lstep.rows - 1) + j))); + float yneg = ((*(c.ptr(Lstep.rows - 2) + j)) + (*(c.ptr(Lstep.rows - 1) + j)))*((*(Ld.ptr(Lstep.rows - 1) + j)) - (*(Ld.ptr(Lstep.rows - 2) + j))); + *(Lstep.ptr(Lstep.rows - 1) + j) = 0.5f*stepsize*(xpos - xneg + ypos - yneg); + } - for (int i = 1; i < Lstep.rows-1; i++) { - float xpos = ((*(c.ptr(i)))+(*(c.ptr(i)+1)))*((*(Ld.ptr(i)+1))-(*(Ld.ptr(i)))); - float xneg = ((*(c.ptr(i)))+(*(c.ptr(i))))*((*(Ld.ptr(i)))-(*(Ld.ptr(i)))); - float ypos = ((*(c.ptr(i)))+(*(c.ptr(i+1))))*((*(Ld.ptr(i+1)))-(*(Ld.ptr(i)))); - float yneg = ((*(c.ptr(i-1)))+(*(c.ptr(i))))*((*(Ld.ptr(i)))-(*(Ld.ptr(i-1)))); - *(Lstep.ptr(i)) = 0.5f*stepsize*(xpos-xneg + ypos-yneg); - } + for (int i = 1; i < Lstep.rows - 1; i++) { + float xpos = ((*(c.ptr(i))) + (*(c.ptr(i)+1)))*((*(Ld.ptr(i)+1)) - (*(Ld.ptr(i)))); + float xneg = ((*(c.ptr(i))) + (*(c.ptr(i))))*((*(Ld.ptr(i))) - (*(Ld.ptr(i)))); + float ypos = ((*(c.ptr(i))) + (*(c.ptr(i + 1))))*((*(Ld.ptr(i + 1))) - (*(Ld.ptr(i)))); + float yneg = ((*(c.ptr(i - 1))) + (*(c.ptr(i))))*((*(Ld.ptr(i))) - (*(Ld.ptr(i - 1)))); + *(Lstep.ptr(i)) = 0.5f*stepsize*(xpos - xneg + ypos - yneg); + } - for (int i = 1; i < Lstep.rows-1; i++) { - float xpos = ((*(c.ptr(i)+Lstep.cols-1))+(*(c.ptr(i)+Lstep.cols-1)))*((*(Ld.ptr(i)+Lstep.cols-1))-(*(Ld.ptr(i)+Lstep.cols-1))); - float xneg = ((*(c.ptr(i)+Lstep.cols-2))+(*(c.ptr(i)+Lstep.cols-1)))*((*(Ld.ptr(i)+Lstep.cols-1))-(*(Ld.ptr(i)+Lstep.cols-2))); - float ypos = ((*(c.ptr(i)+Lstep.cols-1))+(*(c.ptr(i+1)+Lstep.cols-1)))*((*(Ld.ptr(i+1)+Lstep.cols-1))-(*(Ld.ptr(i)+Lstep.cols-1))); - float yneg = ((*(c.ptr(i-1)+Lstep.cols-1))+(*(c.ptr(i)+Lstep.cols-1)))*((*(Ld.ptr(i)+Lstep.cols-1))-(*(Ld.ptr(i-1)+Lstep.cols-1))); - *(Lstep.ptr(i)+Lstep.cols-1) = 0.5f*stepsize*(xpos-xneg + ypos-yneg); - } + for (int i = 1; i < Lstep.rows - 1; i++) { + float xpos = ((*(c.ptr(i)+Lstep.cols - 1)) + (*(c.ptr(i)+Lstep.cols - 1)))*((*(Ld.ptr(i)+Lstep.cols - 1)) - (*(Ld.ptr(i)+Lstep.cols - 1))); + float xneg = ((*(c.ptr(i)+Lstep.cols - 2)) + (*(c.ptr(i)+Lstep.cols - 1)))*((*(Ld.ptr(i)+Lstep.cols - 1)) - (*(Ld.ptr(i)+Lstep.cols - 2))); + float ypos = ((*(c.ptr(i)+Lstep.cols - 1)) + (*(c.ptr(i + 1) + Lstep.cols - 1)))*((*(Ld.ptr(i + 1) + Lstep.cols - 1)) - (*(Ld.ptr(i)+Lstep.cols - 1))); + float yneg = ((*(c.ptr(i - 1) + Lstep.cols - 1)) + (*(c.ptr(i)+Lstep.cols - 1)))*((*(Ld.ptr(i)+Lstep.cols - 1)) - (*(Ld.ptr(i - 1) + Lstep.cols - 1))); + *(Lstep.ptr(i)+Lstep.cols - 1) = 0.5f*stepsize*(xpos - xneg + ypos - yneg); + } - Ld = Ld + Lstep; -} - -//************************************************************************************* -//************************************************************************************* - -/** - * @brief This function checks if a given pixel is a maximum in a local neighbourhood - * @param img Input image where we will perform the maximum search - * @param dsize Half size of the neighbourhood - * @param value Response value at (x,y) position - * @param row Image row coordinate - * @param col Image column coordinate - * @param same_img Flag to indicate if the image value at (x,y) is in the input image - * @return 1->is maximum, 0->otherwise - */ -bool check_maximum_neighbourhood(const cv::Mat& img, int dsize, float value, - int row, int col, bool same_img) { - - bool response = true; - - for (int i = row-dsize; i <= row+dsize; i++) { - for (int j = col-dsize; j <= col+dsize; j++) { - if (i >= 0 && i < img.rows && j >= 0 && j < img.cols) { - if (same_img == true) { - if (i != row || j != col) { - if ((*(img.ptr(i)+j)) > value) { - response = false; - return response; + Ld = Ld + Lstep; } - } - } - else { - if ((*(img.ptr(i)+j)) > value) { - response = false; - return response; - } - } - } - } - } - return response; -} + //************************************************************************************* + //************************************************************************************* + + /** + * @brief This function checks if a given pixel is a maximum in a local neighbourhood + * @param img Input image where we will perform the maximum search + * @param dsize Half size of the neighbourhood + * @param value Response value at (x,y) position + * @param row Image row coordinate + * @param col Image column coordinate + * @param same_img Flag to indicate if the image value at (x,y) is in the input image + * @return 1->is maximum, 0->otherwise + */ + bool check_maximum_neighbourhood(const cv::Mat& img, int dsize, float value, + int row, int col, bool same_img) { + + bool response = true; + + for (int i = row - dsize; i <= row + dsize; i++) { + for (int j = col - dsize; j <= col + dsize; j++) { + if (i >= 0 && i < img.rows && j >= 0 && j < img.cols) { + if (same_img == true) { + if (i != row || j != col) { + if ((*(img.ptr(i)+j)) > value) { + response = false; + return response; + } + } + } + else { + if ((*(img.ptr(i)+j)) > value) { + response = false; + return response; + } + } + } + } + } + + return response; + } + } + } +} \ No newline at end of file diff --git a/modules/features2d/src/kaze/nldiffusion_functions.h b/modules/features2d/src/kaze/nldiffusion_functions.h index d0ece89571..e9d5f03670 100644 --- a/modules/features2d/src/kaze/nldiffusion_functions.h +++ b/modules/features2d/src/kaze/nldiffusion_functions.h @@ -1,4 +1,3 @@ - /** * @file nldiffusion_functions.h * @brief Functions for non-linear diffusion applications: @@ -9,43 +8,40 @@ * @author Pablo F. Alcantarilla */ -#ifndef NLDIFFUSION_FUNCTIONS_H_ -#define NLDIFFUSION_FUNCTIONS_H_ - -//****************************************************************************** -//****************************************************************************** +#ifndef KAZE_NLDIFFUSION_FUNCTIONS_H +#define KAZE_NLDIFFUSION_FUNCTIONS_H // Includes -#include "config.h" +#include "precomp.hpp" //************************************************************************************* //************************************************************************************* -// Gaussian 2D convolution -void gaussian_2D_convolution(const cv::Mat& src, cv::Mat& dst, - int ksize_x, int ksize_y, float sigma); +namespace cv { + namespace details { + namespace kaze { -// Diffusivity functions -void pm_g1(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k); -void pm_g2(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k); -void weickert_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k); -float compute_k_percentile(const cv::Mat& img, float perc, float gscale, - int nbins, int ksize_x, int ksize_y); + // Gaussian 2D convolution + void gaussian_2D_convolution(const cv::Mat& src, cv::Mat& dst, int ksize_x, int ksize_y, float sigma); -// Image derivatives -void compute_scharr_derivatives(const cv::Mat& src, cv::Mat& dst, - int xorder, int yorder, int scale); -void compute_derivative_kernels(cv::OutputArray _kx, cv::OutputArray _ky, - int dx, int dy, int scale); + // Diffusivity functions + void pm_g1(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k); + void pm_g2(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k); + void weickert_diffusivity(const cv::Mat& Lx, const cv::Mat& Ly, cv::Mat& dst, float k); + float compute_k_percentile(const cv::Mat& img, float perc, float gscale, int nbins, int ksize_x, int ksize_y); -// Nonlinear diffusion filtering scalar step -void nld_step_scalar(cv::Mat& Ld, const cv::Mat& c, cv::Mat& Lstep, float stepsize); + // Image derivatives + void compute_scharr_derivatives(const cv::Mat& src, cv::Mat& dst, int xorder, int yorder, int scale); + void compute_derivative_kernels(cv::OutputArray _kx, cv::OutputArray _ky, int dx, int dy, int scale); -// For non-maxima suppresion -bool check_maximum_neighbourhood(const cv::Mat& img, int dsize, float value, - int row, int col, bool same_img); + // Nonlinear diffusion filtering scalar step + void nld_step_scalar(cv::Mat& Ld, const cv::Mat& c, cv::Mat& Lstep, float stepsize); -//************************************************************************************* -//************************************************************************************* + // For non-maxima suppresion + bool check_maximum_neighbourhood(const cv::Mat& img, int dsize, float value, int row, int col, bool same_img); -#endif // NLDIFFUSION_FUNCTIONS_H_ + } + } +} + +#endif