From ad50964f78c43d8af519a1bfd4bafdcba39f8fee Mon Sep 17 00:00:00 2001 From: Anand Mahesh Date: Wed, 6 May 2026 13:44:47 +0530 Subject: [PATCH] Merge pull request #28879 from manand881:feat/akaze-ocl-performance Feat: Add OpenCL support for AKAZE features #28879 Implement OpenCL-accelerated AKAZE feature detection and descriptor extraction Benchmarked on RTX 5060 Ti: 1.2x to 3.31x faster for image sizes from 640x480 to 3840x2160 - Added 8 OpenCL kernels for feature detection and descriptor extraction - Refactored compute_kcontrast to use OpenCV magnitude() for consistency - Added comprehensive tests with >95% accuracy vs CPU baseline ### Pull Request Readiness Checklist See details at https://github.com/opencv/opencv/wiki/How_to_contribute#making-a-good-pull-request - [x] I agree to contribute to the project under Apache 2 License. - [x] To the best of my knowledge, the proposed patch is not based on a code under GPL or another license that is incompatible with OpenCV - [x] The PR is proposed to the proper branch - [ ] There is a reference to the original bug report and related work - [x] There is accuracy test, performance test and test data in opencv_extra repository, if applicable Patch to opencv_extra has the same branch name. - [x] The feature is well documented and sample code can be built with the project CMake --- modules/features2d/perf/opencl/perf_akaze.cpp | 127 +++ modules/features2d/src/akaze.cpp | 42 +- modules/features2d/src/kaze.cpp | 60 ++ modules/features2d/src/kaze/AKAZEFeatures.cpp | 852 ++++++++++++------ modules/features2d/src/kaze/AKAZEFeatures.h | 6 + modules/features2d/src/kaze/KAZEFeatures.cpp | 278 +++++- modules/features2d/src/kaze/KAZEFeatures.h | 10 +- modules/features2d/src/kaze/TEvolution.h | 20 +- .../features2d/src/kaze/akaze_diffusion.hpp | 218 +++++ modules/features2d/src/opencl/akaze.cl | 524 +++++++++++ modules/features2d/src/opencl/kaze.cl | 296 ++++++ modules/features2d/test/ocl/test_akaze.cpp | 140 +++ modules/features2d/test/test_akaze.cpp | 6 +- 13 files changed, 2293 insertions(+), 286 deletions(-) create mode 100644 modules/features2d/perf/opencl/perf_akaze.cpp create mode 100644 modules/features2d/src/kaze/akaze_diffusion.hpp create mode 100644 modules/features2d/src/opencl/kaze.cl create mode 100644 modules/features2d/test/ocl/test_akaze.cpp diff --git a/modules/features2d/perf/opencl/perf_akaze.cpp b/modules/features2d/perf/opencl/perf_akaze.cpp new file mode 100644 index 0000000000..c1b3d1dd9a --- /dev/null +++ b/modules/features2d/perf/opencl/perf_akaze.cpp @@ -0,0 +1,127 @@ +// 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 "../perf_precomp.hpp" +#include "opencv2/ts/ocl_perf.hpp" + +#ifdef HAVE_OPENCL + +namespace opencv_test { +namespace ocl { + +typedef Size_MatType AKAZEFixture; + +OCL_PERF_TEST_P(AKAZEFixture, detectAndCompute, ::testing::Combine(OCL_PERF_ENUM(OCL_SIZE_1, OCL_SIZE_2, OCL_SIZE_3), OCL_PERF_ENUM((MatType)CV_8UC1))) +{ + const Size_MatType_t params = GetParam(); + const Size srcSize = get<0>(params); + const int type = get<1>(params); + + checkDeviceMaxMemoryAllocSize(srcSize, type); + + UMat img(srcSize, type), mask; + declare.in(img, WARMUP_RNG); + + Ptr akaze = AKAZE::create(AKAZE::DESCRIPTOR_MLDB_UPRIGHT, 0, 3, 0.001f, 1, 1, KAZE::DIFF_PM_G2); + vector points; + UMat descriptors; + + OCL_TEST_CYCLE() akaze->detectAndCompute(img, mask, points, descriptors, false); + + EXPECT_GT(points.size(), 20u); + EXPECT_EQ((size_t)descriptors.rows, points.size()); + SANITY_CHECK_NOTHING(); +} + +typedef Size_MatType AKAZEOclMldbUprightFixture; + +OCL_PERF_TEST_P(AKAZEOclMldbUprightFixture, detectAndComputeMLDB, + ::testing::Combine( + ::testing::Values( + cv::Size(320, 240), cv::Size(640, 480), cv::Size(960, 540), + cv::Size(1280, 720), cv::Size(1920, 1080), cv::Size(2560, 1440), cv::Size(3840, 2160) + ), + OCL_PERF_ENUM((MatType)CV_8UC1) + )) +{ + const Size_MatType_t params = GetParam(); + const cv::Size srcSize = get<0>(params); + + UMat img(srcSize, CV_8U); + declare.in(img, WARMUP_RNG); + + Ptr det = AKAZE::create(AKAZE::DESCRIPTOR_MLDB_UPRIGHT, 0, 3, 0.001f, 4, 4, KAZE::DIFF_PM_G2); + vector kps; + UMat descriptors; + + OCL_TEST_CYCLE() det->detectAndCompute(img, noArray(), kps, descriptors, false); + + EXPECT_GT(kps.size(), 0u); + SANITY_CHECK_NOTHING(); +} + +typedef tuple AKAZEOclRealImagesParams; +typedef TestBaseWithParam AKAZEOclRealImagesFixture; + +OCL_PERF_TEST_P(AKAZEOclRealImagesFixture, detectAndComputeMLDBRealImages, + ::testing::Combine( + ::testing::Values( + std::string("stitching/boat1.jpg"), + std::string("stitching/boat2.jpg"), + std::string("stitching/boat3.jpg"), + std::string("stitching/boat4.jpg"), + std::string("stitching/boat5.jpg"), + std::string("stitching/boat6.jpg") + ), + ::testing::Values(0) + )) +{ + const std::string imgName = get<0>(GetParam()); + + Mat img_mat = imread(getDataPath(imgName), IMREAD_GRAYSCALE); + if (img_mat.empty()) + throw cvtest::SkipTestException("Image not found: " + imgName); + + UMat img; + img_mat.copyTo(img); + declare.in(img); + + Ptr det = AKAZE::create(AKAZE::DESCRIPTOR_MLDB_UPRIGHT, 0, 3, 0.001f, 4, 4, KAZE::DIFF_PM_G2); + vector kps; + UMat descriptors; + + OCL_TEST_CYCLE() det->detectAndCompute(img, noArray(), kps, descriptors, false); + + EXPECT_GT(kps.size(), 0u); + SANITY_CHECK_NOTHING(); +} + +typedef Size_MatType AKAZEOclKazeUprightFixture; + +OCL_PERF_TEST_P(AKAZEOclKazeUprightFixture, detectAndComputeKAZEUpright, + ::testing::Combine( + ::testing::Values(cv::Size(320, 240), cv::Size(640, 480), cv::Size(1280, 720)), + OCL_PERF_ENUM((MatType)CV_8UC1) + )) +{ + const Size_MatType_t params = GetParam(); + const cv::Size srcSize = get<0>(params); + + UMat img(srcSize, CV_8U); + declare.in(img, WARMUP_RNG); + + Ptr det = AKAZE::create(AKAZE::DESCRIPTOR_KAZE_UPRIGHT, 0, 3, 0.001f, 4, 4, KAZE::DIFF_PM_G2); + vector kps; + UMat descriptors; + + OCL_TEST_CYCLE() det->detectAndCompute(img, noArray(), kps, descriptors, false); + + EXPECT_GT(kps.size(), 0u); + SANITY_CHECK_NOTHING(); +} + +} // ocl +} // opencv_test + +#endif // HAVE_OPENCL diff --git a/modules/features2d/src/akaze.cpp b/modules/features2d/src/akaze.cpp index a41ee55200..d7868959d1 100644 --- a/modules/features2d/src/akaze.cpp +++ b/modules/features2d/src/akaze.cpp @@ -187,11 +187,36 @@ namespace cv options.diffusivity = diffusivity; AKAZEFeatures impl(options); - impl.Create_Nonlinear_Scale_Space(image); - if (!useProvidedKeypoints) + UMatPyramid uPyr; +#ifdef HAVE_OPENCL + // NOTE: The AKAZE kernels use excessive private memory 336 bytes per work-item: + // float histogram[36] and float values[48] arrays. On 32-bit Windows, this + // exceeds the driver's private memory limits, causing SEH exception (0xc0000005). +#if defined(_M_IX86) || defined(__i386__) + bool use_opencl = false; +#else + bool use_opencl = cv::ocl::useOpenCL() && image.isUMat() && descriptor == DESCRIPTOR_MLDB_UPRIGHT + && !ocl::Device::getDefault().hostUnifiedMemory(); +#endif + if (use_opencl) { - impl.Feature_Detection(keypoints); + impl.GetEvolutionPyramid(uPyr); // Get initialized pyramid + impl.Create_Nonlinear_Scale_Space_UMat(image, uPyr); + if (!useProvidedKeypoints) + { + impl.Feature_Detection_UMat(uPyr, keypoints); + } + } + else +#endif + { + impl.Create_Nonlinear_Scale_Space(image); + + if (!useProvidedKeypoints) + { + impl.Feature_Detection(keypoints); + } } if (!mask.empty()) @@ -207,7 +232,16 @@ namespace cv if(descriptors.needed()) { - impl.Compute_Descriptors(keypoints, descriptors); +#ifdef HAVE_OPENCL + if (use_opencl) + { + impl.Compute_Descriptors_UMat(keypoints, descriptors, uPyr); + } + else +#endif + { + impl.Compute_Descriptors(keypoints, descriptors); + } CV_Assert((descriptors.empty() || descriptors.cols() == descriptorSize())); CV_Assert((descriptors.empty() || (descriptors.type() == descriptorType()))); diff --git a/modules/features2d/src/kaze.cpp b/modules/features2d/src/kaze.cpp index 881c5aab1e..30cb304ab4 100644 --- a/modules/features2d/src/kaze.cpp +++ b/modules/features2d/src/kaze.cpp @@ -112,6 +112,66 @@ namespace cv { CV_INSTRUMENT_REGION(); +#ifdef HAVE_OPENCL + bool use_opencl = ocl::useOpenCL() && upright + && !ocl::Device::getDefault().hostUnifiedMemory(); + if (use_opencl) + { + UMat img_gray = image.getUMat(); + if (img_gray.channels() > 1) + { + UMat tmp; + cvtColor(img_gray, tmp, COLOR_BGR2GRAY); + img_gray = tmp; + } + + UMat img1_32; + int depth = img_gray.depth(); + if ( depth == CV_32F ) + img1_32 = img_gray; + else if ( depth == CV_8U ) + img_gray.convertTo(img1_32, CV_32F, 1.0 / 255.0, 0); + else if ( depth == CV_16U ) + img_gray.convertTo(img1_32, CV_32F, 1.0 / 65535.0, 0); + + CV_Assert( ! img1_32.empty() ); + + KAZEOptions options; + options.img_width = img1_32.cols; + options.img_height = img1_32.rows; + options.extended = extended; + options.upright = upright; + options.dthreshold = threshold; + options.omax = octaves; + options.nsublevels = sublevels; + options.diffusivity = diffusivity; + + KAZEFeatures impl(options); + UTPyramid uPyr; + impl.Create_Nonlinear_Scale_Space_UMat(img1_32, uPyr); + + if (!useProvidedKeypoints) + { + impl.Feature_Detection_UMat(uPyr, keypoints); + } + + if (!mask.empty()) + { + cv::KeyPointsFilter::runByPixelsMask(keypoints, mask.getMat()); + } + + if( descriptors.needed() ) + { + impl.Compute_Descriptors_UMat(keypoints, descriptors, uPyr); + + CV_Assert((!descriptors.empty() || descriptors.cols() == descriptorSize())); + CV_Assert((!descriptors.empty() || (descriptors.type() == descriptorType()))); + } + return; + } +#endif + + // CPU path cv::Mat img = image.getMat(); if (img.channels() > 1) cvtColor(image, img, COLOR_BGR2GRAY); diff --git a/modules/features2d/src/kaze/AKAZEFeatures.cpp b/modules/features2d/src/kaze/AKAZEFeatures.cpp index e91bb5e349..83cf6185be 100644 --- a/modules/features2d/src/kaze/AKAZEFeatures.cpp +++ b/modules/features2d/src/kaze/AKAZEFeatures.cpp @@ -12,6 +12,7 @@ #include "nldiffusion_functions.h" #include "utils.h" #include "opencl_kernels_features2d.hpp" +#include "akaze_diffusion.hpp" #include @@ -111,185 +112,6 @@ static inline int getGaussianKernelSize(float sigma) { 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(0) + 1; /* Skip the left-most column by +1 */ - lf_c = Lf.ptr(0) + 1; - lt_b = Lt.ptr(1) + 1; - lf_b = Lf.ptr(1) + 1; - - // fill the corner to prevent uninitialized values - dst = Lstep.ptr(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(row - 1); - lf_a = Lf.ptr(row - 1); - lt_c = Lt.ptr(row ); - lf_c = Lf.ptr(row ); - lt_b = Lt.ptr(row + 1); - lf_b = Lf.ptr(row + 1); - dst = Lstep.ptr(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(row - 1) + 1; /* Skip the left-most column by +1 */ - lf_a = Lf.ptr(row - 1) + 1; - lt_c = Lt.ptr(row ) + 1; - lf_c = Lf.ptr(row ) + 1; - - // fill the corner to prevent uninitialized values - dst = Lstep.ptr(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 CV_OVERRIDE - { - 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(InputArray Lt_, InputArray Lf_, OutputArray Lstep_, float step_size) -{ - if(!Lt_.isContinuous()) - return false; - - UMat Lt = Lt_.getUMat(); - UMat Lf = Lf_.getUMat(); - UMat Lstep = Lstep_.getUMat(); - - 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(InputArray Lt_, InputArray Lf_, OutputArray Lstep_, float step_size) -{ - CV_INSTRUMENT_REGION(); - - Lstep_.create(Lt_.size(), Lt_.type()); - - CV_OCL_RUN(Lt_.isUMat() && Lf_.isUMat() && Lstep_.isUMat(), - ocl_non_linear_diffusion_step(Lt_, Lf_, Lstep_, step_size)); - - Mat Lt = Lt_.getMat(); - Mat Lf = Lf_.getMat(); - Mat Lstep = Lstep_.getMat(); - parallel_for_(Range(0, Lt.rows), NonLinearScalarDiffusionStep(Lt, Lf, Lstep, 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 @@ -310,38 +132,42 @@ compute_kcontrast(InputArray Lx_, InputArray Ly_, float perc, int nbins) Mat Lx = Lx_.getMat(); Mat Ly = Ly_.getMat(); - // 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(); + int rows = Lx.rows - 2; + int cols = Lx.cols - 2; + const int total = rows * cols; + + if (total <= 0) + return 0.03f; + + Mat Lx_inner = Lx(Rect(1, 1, cols, rows)); + Mat Ly_inner = Ly(Rect(1, 1, cols, rows)); + Mat modgs; + magnitude(Lx_inner, Ly_inner, modgs); + float hmax = 0.0f; - - for (int i = 1; i < Lx.rows - 1; i++) { - const float *lx = Lx.ptr(i) + 1; - const float *ly = Ly.ptr(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); - } + const float* modg_ptr = modgs.ptr(); + for (int i = 0; i < total; i++) { + if (modg_ptr[i] > hmax) hmax = modg_ptr[i]; } - modg = modgs.ptr(); if (hmax == 0.0f) - return 0.03f; // e.g. a blank image + return 0.03f; - // Compute the bin numbers: the value range [0, hmax] -> [0, nbins-1] - modgs *= (nbins - 1) / hmax; + Mat modgs_mutable = modgs.clone(); + float* modg_ptr_w = modgs_mutable.ptr(); + float scale = (nbins - 1) / hmax; + for (int i = 0; i < total; i++) { + modg_ptr_w[i] = modg_ptr[i] * scale; + } - // Count up histogram std::vector hist(nbins, 0); - for (int i = 0; i < total; i++) - hist[(int)modg[i]]++; + for (int i = 0; i < total; i++) { + int bin = (int)modg_ptr_w[i]; + if (bin >= 0 && bin < nbins) + hist[bin]++; + } - // Now find the perc of the histogram percentile - const int nthreshold = (int)((total - hist[0]) * perc); // Exclude hist[0] as background + const int nthreshold = (int)((total - hist[0]) * perc); int nelements = 0; for (int k = 1; k < nbins; k++) { if (nelements >= nthreshold) @@ -355,72 +181,102 @@ compute_kcontrast(InputArray Lx_, InputArray Ly_, float perc, int nbins) #ifdef HAVE_OPENCL static inline bool -ocl_pm_g2(InputArray Lx_, InputArray Ly_, OutputArray Lflow_, float kcontrast) +ocl_compute_kcontrast(InputArray Lx_, InputArray Ly_, float perc, int nbins, float& kcontrast_out) { UMat Lx = Lx_.getUMat(); UMat Ly = Ly_.getUMat(); - UMat Lflow = Lflow_.getUMat(); - int total = Lx.rows * Lx.cols; - size_t globalSize[] = {(size_t)total}; + int rows = Lx.rows - 2; + int cols = Lx.cols - 2; + int total = rows * cols; - ocl::Kernel ker("AKAZE_pm_g2", ocl::features2d::akaze_oclsrc); - if( ker.empty() ) + if (total <= 0) return false; - return ker.args( - ocl::KernelArg::PtrReadOnly(Lx), - ocl::KernelArg::PtrReadOnly(Ly), - ocl::KernelArg::PtrWriteOnly(Lflow), - kcontrast, total).run(1, globalSize, 0, true); + UMat modgs(rows, cols, CV_32F); + UMat Lx_inner = Lx(Rect(1, 1, cols, rows)); + UMat Ly_inner = Ly(Rect(1, 1, cols, rows)); + + magnitude(Lx_inner, Ly_inner, modgs); + + // Download once and compute histogram on CPU + Mat modgs_mat; + modgs.copyTo(modgs_mat); + + const float* ptr = modgs_mat.ptr(); + float hmax = 0.0f; + for (int i = 0; i < total; i++) { + if (ptr[i] > hmax) hmax = ptr[i]; + } + + if (hmax == 0.0f) { + kcontrast_out = 0.03f; + return true; + } + + float scale = (nbins - 1) / hmax; + float* ptr_w = (float*)ptr; // We already have a writable copy via copyTo + for (int i = 0; i < total; i++) { + ptr_w[i] = ptr[i] * scale; + } + + Mat hist_mat(1, nbins, CV_32SC1, Scalar(0)); + for (int i = 0; i < total; i++) { + int bin = (int)ptr_w[i]; + if (bin >= 0 && bin < nbins) { + hist_mat.at(0, bin)++; + } + } + + const int nthreshold = (int)((total - hist_mat.at(0, 0)) * perc); + int nelements = 0; + for (int k = 1; k < nbins; k++) { + if (nelements >= nthreshold) { + kcontrast_out = (float)hmax * k / nbins; + return true; + } + nelements += hist_mat.at(0, k); + } + + kcontrast_out = 0.03f; + return true; } #endif // HAVE_OPENCL -static inline void -compute_diffusivity(InputArray Lx, InputArray Ly, OutputArray Lflow, float kcontrast, KAZE::DiffusivityType 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(Lx.isUMat() && Ly.isUMat() && Lflow.isUMat(), 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_(Error::StsError, ("Diffusivity is not supported: %d", static_cast(diffusivity))); - break; - } -} - /** * @brief Converts input image to grayscale float image * - * @param image any image - * @param dst grayscale float image + * @param image any image (Mat or UMat) + * @param dst grayscale float image (Mat or UMat) */ static inline void prepareInputImage(InputArray image, OutputArray dst) { - Mat img = image.getMat(); - if (img.channels() > 1) - cvtColor(image, img, COLOR_BGR2GRAY); + if (image.isUMat()) { + UMat img = image.getUMat(); + if (img.channels() > 1) { + UMat tmp; + cvtColor(image, tmp, COLOR_BGR2GRAY); + img = tmp; + } - if ( img.depth() == CV_32F ) - dst.assign(img); - else if ( img.depth() == CV_8U ) - img.convertTo(dst, CV_32F, 1.0 / 255.0, 0); - else if ( img.depth() == CV_16U ) - img.convertTo(dst, CV_32F, 1.0 / 65535.0, 0); + if (img.depth() == CV_32F) + dst.assign(img); + else if (img.depth() == CV_8U) + img.convertTo(dst, CV_32F, 1.0 / 255.0, 0); + else if (img.depth() == CV_16U) + img.convertTo(dst, CV_32F, 1.0 / 65535.0, 0); + } else { + Mat img = image.getMat(); + if (img.channels() > 1) + cvtColor(image, img, COLOR_BGR2GRAY); + + if (img.depth() == CV_32F) + dst.assign(img); + else if (img.depth() == CV_8U) + img.convertTo(dst, CV_32F, 1.0 / 255.0, 0); + else if (img.depth() == CV_16U) + img.convertTo(dst, CV_32F, 1.0 / 65535.0, 0); + } } /** @@ -458,8 +314,17 @@ create_nonlinear_scale_space(InputArray image, const AKAZEOptions &options, 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, Ly, options.kcontrast_percentile, options.kcontrast_nbins); + + // compute the kcontrast factor - try OCL path first if working with UMat + float kcontrast = 0.0f; +#ifdef HAVE_OPENCL + bool kcontrast_ocl_success = false; + if (std::is_same::value) { + kcontrast_ocl_success = ocl_compute_kcontrast(Lx, Ly, options.kcontrast_percentile, options.kcontrast_nbins, kcontrast); + } + if (!kcontrast_ocl_success) +#endif + kcontrast = compute_kcontrast(Lx, Ly, options.kcontrast_percentile, options.kcontrast_nbins); // Now generate the rest of evolution levels for (size_t i = 1; i < evolution.size(); i++) { @@ -520,7 +385,7 @@ convertScalePyramid(const std::vector >& src, std::vector< */ void AKAZEFeatures::Create_Nonlinear_Scale_Space(InputArray image) { - if (ocl::isOpenCLActivated() && image.isUMat()) { + if (ocl::isOpenCLActivated() && image.isUMat() && !ocl::Device::getDefault().hostUnifiedMemory()) { // will run OCL version of scale space pyramid UMatPyramid uPyr; // init UMat pyramid with sizes @@ -534,6 +399,27 @@ void AKAZEFeatures::Create_Nonlinear_Scale_Space(InputArray image) } } +#ifdef HAVE_OPENCL +/** + * @brief Get initialized evolution pyramid for UMatPyramid + * @param uPyr Output UMatPyramid initialized with sizes from evolution_ + */ +void AKAZEFeatures::GetEvolutionPyramid(UMatPyramid& uPyr) +{ + convertScalePyramid(evolution_, uPyr); +} + +/** + * @brief This method creates the nonlinear scale space for a given image (UMatPyramid version) + * @param image Input image for which the nonlinear scale space needs to be created + * @param uPyr UMatPyramid to store the evolution + */ +void AKAZEFeatures::Create_Nonlinear_Scale_Space_UMat(InputArray image, UMatPyramid& uPyr) +{ + create_nonlinear_scale_space(image, options_, tsteps_, uPyr); +} +#endif // HAVE_OPENCL + /* ************************************************************************* */ #ifdef HAVE_OPENCL @@ -558,7 +444,7 @@ ocl_compute_determinant(InputArray Lxx_, InputArray Lxy_, InputArray Lyy_, ocl::KernelArg::PtrReadOnly(Lxy), ocl::KernelArg::PtrReadOnly(Lyy), ocl::KernelArg::PtrWriteOnly(Ldet), - sigma, total).run(1, globalSize, 0, true); + sigma, total).run(1, globalSize, 0, false); } #endif // HAVE_OPENCL @@ -714,6 +600,202 @@ find_neighbor_point(const int x, const int y, const Mat &mask, const int search_ return false; } +#ifdef HAVE_OPENCL +static inline bool +ocl_find_extrema_same_scale(InputArray Ldet_, OutputArray keypoint_mask_, float threshold, int border) +{ + UMat Ldet = Ldet_.getUMat(); + UMat keypoint_mask = keypoint_mask_.getUMat(); + + int rows = Ldet.rows; + int cols = Ldet.cols; + + size_t globalSize[] = {(size_t)cols, (size_t)rows}; + + ocl::Kernel ker("AKAZE_find_extrema_same_scale", ocl::features2d::akaze_oclsrc); + if( ker.empty() ) + return false; + + bool success = ker.set(0, ocl::KernelArg::ReadOnly(Ldet)); + success = success && ker.set(1, rows); + success = success && ker.set(2, cols); + success = success && ker.set(3, threshold); + success = success && ker.set(4, border); + success = success && ker.set(5, ocl::KernelArg::PtrWriteOnly(keypoint_mask)); + + if (!success) + return false; + + return ker.run(2, globalSize, 0, false); +} + +static inline bool +ocl_cross_scale_filter_lower(InputArray keypoints_current_, InputOutputArray keypoints_lower_, + InputArray ldet_current_, InputArray ldet_lower_, int diff_ratio, int search_radius) +{ + UMat keypoints_current = keypoints_current_.getUMat(); + UMat keypoints_lower = keypoints_lower_.getUMat(); + UMat ldet_current = ldet_current_.getUMat(); + UMat ldet_lower = ldet_lower_.getUMat(); + + int rows = keypoints_current.rows; + int cols = keypoints_current.cols; + int lower_rows = keypoints_lower.rows; + int lower_cols = keypoints_lower.cols; + + size_t globalSize[] = {(size_t)cols, (size_t)rows}; + + ocl::Kernel ker("AKAZE_cross_scale_filter_lower", ocl::features2d::akaze_oclsrc); + if( ker.empty() ) + return false; + + bool success = ker.set(0, ocl::KernelArg::ReadOnly(keypoints_current)); + success = success && ker.set(1, ocl::KernelArg::ReadWrite(keypoints_lower)); + success = success && ker.set(2, ocl::KernelArg::ReadOnly(ldet_current)); + success = success && ker.set(3, ocl::KernelArg::ReadOnly(ldet_lower)); + success = success && ker.set(4, rows); + success = success && ker.set(5, cols); + success = success && ker.set(6, lower_rows); + success = success && ker.set(7, lower_cols); + success = success && ker.set(8, diff_ratio); + success = success && ker.set(9, search_radius); + return success && ker.run(2, globalSize, 0, false); +} + +static inline bool +ocl_cross_scale_filter_upper(InputArray keypoints_current_, InputOutputArray keypoints_upper_, + InputArray ldet_current_, InputArray ldet_upper_, int diff_ratio, int search_radius) +{ + UMat keypoints_current = keypoints_current_.getUMat(); + UMat keypoints_upper = keypoints_upper_.getUMat(); + UMat ldet_current = ldet_current_.getUMat(); + UMat ldet_upper = ldet_upper_.getUMat(); + + int rows = keypoints_current.rows; + int cols = keypoints_current.cols; + int upper_rows = keypoints_upper.rows; + int upper_cols = keypoints_upper.cols; + + size_t globalSize[] = {(size_t)cols, (size_t)rows}; + + ocl::Kernel ker("AKAZE_cross_scale_filter_upper", ocl::features2d::akaze_oclsrc); + if( ker.empty() ) + return false; + + bool success = ker.set(0, ocl::KernelArg::ReadOnly(keypoints_current)); + success = success && ker.set(1, ocl::KernelArg::ReadWrite(keypoints_upper)); + success = success && ker.set(2, ocl::KernelArg::ReadOnly(ldet_current)); + success = success && ker.set(3, ocl::KernelArg::ReadOnly(ldet_upper)); + success = success && ker.set(4, rows); + success = success && ker.set(5, cols); + success = success && ker.set(6, upper_rows); + success = success && ker.set(7, upper_cols); + success = success && ker.set(8, diff_ratio); + success = success && ker.set(9, search_radius); + return success && ker.run(2, globalSize, 0, false); +} + +static inline bool +ocl_subpixel_refinement_orientation(InputArray keypoints_, InputArray ldet_, + InputArray Lx_, InputArray Ly_, + int rows, int cols, + float octave_ratio, float esigma, int octave, int level, + OutputArray output_count_, OutputArray output_x_, OutputArray output_y_, + OutputArray output_size_, OutputArray output_response_, + OutputArray output_octave_, OutputArray output_class_id_, + OutputArray output_angle_, + int max_output) +{ + UMat keypoints = keypoints_.getUMat(); + UMat ldet = ldet_.getUMat(); + UMat Lx = Lx_.getUMat(); + UMat Ly = Ly_.getUMat(); + UMat output_count = output_count_.getUMat(); + UMat output_x = output_x_.getUMat(); + UMat output_y = output_y_.getUMat(); + UMat output_size = output_size_.getUMat(); + UMat output_response = output_response_.getUMat(); + UMat output_octave = output_octave_.getUMat(); + UMat output_class_id = output_class_id_.getUMat(); + UMat output_angle = output_angle_.getUMat(); + + size_t globalSize[] = {(size_t)cols, (size_t)rows}; + + ocl::Kernel ker("AKAZE_subpixel_refinement_orientation", ocl::features2d::akaze_oclsrc); + if( ker.empty() ) + return false; + + bool success = ker.set(0, ocl::KernelArg::ReadOnly(keypoints)); + success = success && ker.set(1, ocl::KernelArg::ReadOnly(ldet)); + success = success && ker.set(2, ocl::KernelArg::ReadOnly(Lx)); + success = success && ker.set(3, ocl::KernelArg::ReadOnly(Ly)); + success = success && ker.set(4, rows); + success = success && ker.set(5, cols); + success = success && ker.set(6, octave_ratio); + success = success && ker.set(7, esigma); + success = success && ker.set(8, octave); + success = success && ker.set(9, level); + success = success && ker.set(10, ocl::KernelArg::PtrWriteOnly(output_count)); + success = success && ker.set(11, ocl::KernelArg::PtrWriteOnly(output_x)); + success = success && ker.set(12, ocl::KernelArg::PtrWriteOnly(output_y)); + success = success && ker.set(13, ocl::KernelArg::PtrWriteOnly(output_size)); + success = success && ker.set(14, ocl::KernelArg::PtrWriteOnly(output_response)); + success = success && ker.set(15, ocl::KernelArg::PtrWriteOnly(output_octave)); + success = success && ker.set(16, ocl::KernelArg::PtrWriteOnly(output_class_id)); + success = success && ker.set(17, ocl::KernelArg::PtrWriteOnly(output_angle)); + success = success && ker.set(18, max_output); + + if (!success) + return false; + + return ker.run(2, globalSize, 0, false); +} + +static inline bool +ocl_compute_mldb_descriptor_level(InputArray keypoints_x_, InputArray keypoints_y_, + InputArray keypoints_size_, InputArray keypoints_angle_, + InputArray Lx_, InputArray Ly_, InputArray Lt_, + int rows, int cols, float octave_ratio, + int descriptor_size, OutputArray output_descriptors_) +{ + UMat keypoints_x = keypoints_x_.getUMat(); + UMat keypoints_y = keypoints_y_.getUMat(); + UMat keypoints_size = keypoints_size_.getUMat(); + UMat keypoints_angle = keypoints_angle_.getUMat(); + UMat Lx = Lx_.getUMat(); + UMat Ly = Ly_.getUMat(); + UMat Lt = Lt_.getUMat(); + UMat output_descriptors = output_descriptors_.getUMat(); + + int num_keypoints = (int)keypoints_x.total(); + + size_t globalSize[] = {(size_t)num_keypoints}; + + ocl::Kernel ker("AKAZE_compute_mldb_descriptor_level", ocl::features2d::akaze_oclsrc); + if( ker.empty() ) + return false; + + bool success = ker.set(0, ocl::KernelArg::ReadOnly(keypoints_x)); + success = success && ker.set(1, ocl::KernelArg::ReadOnly(keypoints_y)); + success = success && ker.set(2, ocl::KernelArg::ReadOnly(keypoints_size)); + success = success && ker.set(3, ocl::KernelArg::ReadOnly(keypoints_angle)); + success = success && ker.set(4, num_keypoints); + success = success && ker.set(5, ocl::KernelArg::ReadOnly(Lx)); + success = success && ker.set(6, ocl::KernelArg::ReadOnly(Ly)); + success = success && ker.set(7, ocl::KernelArg::ReadOnly(Lt)); + success = success && ker.set(8, rows); + success = success && ker.set(9, cols); + success = success && ker.set(10, octave_ratio); + success = success && ker.set(11, descriptor_size); + success = success && ker.set(12, ocl::KernelArg::PtrWriteOnly(output_descriptors)); + + if (!success) + return false; + + return ker.run(1, globalSize, 0, false); +} +#endif // HAVE_OPENCL + /** * @brief Find keypoints in parallel for each pyramid layer */ @@ -785,6 +867,264 @@ private: float dthreshold_; ///< Detector response threshold to accept point }; +#ifdef HAVE_OPENCL +/** + * @brief Find keypoints in parallel for each pyramid layer (OpenCL version) + */ +class FindKeypointsSameScaleOCL : public ParallelLoopBody +{ +public: + explicit FindKeypointsSameScaleOCL(const UMatPyramid& ev, + std::vector& kpts, float dthreshold) + : evolution_(&ev), keypoints_by_layers_(&kpts), dthreshold_(dthreshold) + {} + + void operator()(const Range& range) const CV_OVERRIDE + { + for (int i = range.start; i < range.end; i++) + { + const Evolution &e = (*evolution_)[i]; + UMat &kpts = (*keypoints_by_layers_)[i]; + // this mask will hold positions of keypoints in this level + kpts = UMat::zeros(e.Ldet.size(), CV_8UC1); + + // if border is too big we shouldn't search any keypoints + if (e.border + 1 >= e.Ldet.rows) + continue; + + // Use OpenCL kernel for same-scale extrema detection only + ocl_find_extrema_same_scale(e.Ldet, kpts, dthreshold_, e.border); + } + } + +private: + const UMatPyramid* evolution_; + std::vector* keypoints_by_layers_; + float dthreshold_; +}; +#endif // HAVE_OPENCL + +#ifdef HAVE_OPENCL +/** + * @brief This method finds extrema in the nonlinear scale space (UMatPyramid version) + * @param keypoints_by_layers Output vectors of detected keypoints; one vector for each evolution level + */ +static inline void +Find_Scale_Space_Extrema_UMat(UMatPyramid& evolution, std::vector& keypoints_by_layers, float dthreshold) +{ + CV_INSTRUMENT_REGION(); + + keypoints_by_layers.resize(evolution.size()); + + // find points in the same level using OpenCL + FindKeypointsSameScaleOCL(evolution, keypoints_by_layers, dthreshold) + (Range(0, (int)evolution.size())); +} + +/** + * @brief This method selects interesting keypoints through the nonlinear scale space (OpenCL version) + * @param uPyr UMatPyramid evolution + * @param kpts Vector of detected keypoints + */ +void AKAZEFeatures::Feature_Detection_UMat(UMatPyramid& uPyr, std::vector& kpts) +{ + CV_INSTRUMENT_REGION(); + + kpts.clear(); + std::vector keypoints_by_layers; + Find_Scale_Space_Extrema_UMat(uPyr, keypoints_by_layers, options_.dthreshold); + + // Cross-scale filtering using OpenCL (GPU-only) + for (size_t i = 1; i < keypoints_by_layers.size(); i++) { + const int diff_ratio = (int)uPyr[i].octave_ratio / (int)uPyr[i-1].octave_ratio; + const int search_radius = uPyr[i].sigma_size * diff_ratio; + ocl_cross_scale_filter_lower(keypoints_by_layers[i], keypoints_by_layers[i-1], + uPyr[i].Ldet, uPyr[i-1].Ldet, diff_ratio, search_radius); + } + + for (int i = (int)keypoints_by_layers.size() - 2; i >= 0; i--) { + const int diff_ratio = (int)uPyr[i+1].octave_ratio / (int)uPyr[i].octave_ratio; + const int search_radius = uPyr[i+1].sigma_size; + ocl_cross_scale_filter_upper(keypoints_by_layers[i], keypoints_by_layers[i+1], + uPyr[i].Ldet, uPyr[i+1].Ldet, diff_ratio, search_radius); + } + + // Submit ALL subpixel+orientation kernels async before any download. + // OpenCL's in-order queue guarantees correct ordering relative to the + // cross-scale filter kernels already queued above. + int num_levels = (int)keypoints_by_layers.size(); + std::vector max_kpts_per_level(num_levels, 0); + std::vector level_active(num_levels, false); + std::vector out_count(num_levels), out_x(num_levels), out_y(num_levels); + std::vector out_size(num_levels), out_response(num_levels); + std::vector out_octave(num_levels), out_class_id(num_levels), out_angle(num_levels); + + for (int i = 0; i < num_levels; i++) { + const Evolution& e = uPyr[i]; + if (e.border + 1 >= e.Ldet.rows) + continue; + level_active[i] = true; + int max_kpts = (int)(e.Ldet.total() * 0.1f); + max_kpts_per_level[i] = max_kpts; + out_count[i] = UMat(1, 1, CV_32S, Scalar::all(0)); + out_x[i] = UMat(max_kpts, 1, CV_32F); + out_y[i] = UMat(max_kpts, 1, CV_32F); + out_size[i] = UMat(max_kpts, 1, CV_32F); + out_response[i] = UMat(max_kpts, 1, CV_32F); + out_octave[i] = UMat(max_kpts, 1, CV_32S); + out_class_id[i] = UMat(max_kpts, 1, CV_32S); + out_angle[i] = UMat(max_kpts, 1, CV_32F); + ocl_subpixel_refinement_orientation(keypoints_by_layers[i], e.Ldet, e.Lx, e.Ly, + e.Ldet.rows, e.Ldet.cols, + e.octave_ratio, e.esigma, e.octave, i, + out_count[i], out_x[i], out_y[i], out_size[i], out_response[i], + out_octave[i], out_class_id[i], out_angle[i], max_kpts); + } + + // Single sync point: the first getMat() flushes the entire queue. + // All subsequent getMat() calls on already-done buffers return immediately. + std::vector all_kpts; + for (int i = 0; i < num_levels; i++) { + if (!level_active[i]) continue; + int max_kpts = max_kpts_per_level[i]; + + Mat count_mat = out_count[i].getMat(ACCESS_READ); + int num_kpts_level = std::min(count_mat.at(0, 0), max_kpts); + if (num_kpts_level == 0) continue; + + Mat x_mat = out_x[i].getMat(ACCESS_READ); + Mat y_mat = out_y[i].getMat(ACCESS_READ); + Mat size_mat = out_size[i].getMat(ACCESS_READ); + Mat resp_mat = out_response[i].getMat(ACCESS_READ); + Mat oct_mat = out_octave[i].getMat(ACCESS_READ); + Mat cid_mat = out_class_id[i].getMat(ACCESS_READ); + Mat ang_mat = out_angle[i].getMat(ACCESS_READ); + + int start = (int)all_kpts.size(); + all_kpts.resize(start + num_kpts_level); + for (int j = 0; j < num_kpts_level; j++) { + KeyPoint& kp = all_kpts[start + j]; + kp.pt.x = x_mat.at(j, 0); + kp.pt.y = y_mat.at(j, 0); + kp.size = size_mat.at(j, 0); + kp.response = resp_mat.at(j, 0); + kp.octave = oct_mat.at(j, 0); + kp.class_id = cid_mat.at(j, 0); + kp.angle = ang_mat.at(j, 0); + } + } + + kpts = all_kpts; +} +#endif // HAVE_OPENCL + +#ifdef HAVE_OPENCL +/** + * @brief Compute descriptors using OpenCL (GPU-only) + * @param kpts Input keypoints + * @param desc Output descriptors + * @param uPyr UMatPyramid evolution + */ +void AKAZEFeatures::Compute_Descriptors_UMat(std::vector& kpts, OutputArray desc, UMatPyramid& uPyr) +{ + CV_INSTRUMENT_REGION(); + + if (kpts.empty()) { + desc.create(0, 0, CV_8UC1); + return; + } + + // Determine descriptor size based on options + int descriptor_size = 64; + int descriptor_type = CV_32FC1; + if (options_.descriptor >= AKAZE::DESCRIPTOR_MLDB_UPRIGHT) + { + int descriptor_bits = (options_.descriptor_size == 0) + ? (6 + 36 + 120)*options_.descriptor_channels + : options_.descriptor_size; + descriptor_size = divUp(descriptor_bits, 8); + descriptor_type = CV_8UC1; + } + + desc.create((int)kpts.size(), descriptor_size, descriptor_type); + UMat descriptors = desc.getUMat(); + + // For MLDB_UPRIGHT descriptor, use GPU computation (MLDB with rotation falls back to CPU) + if (options_.descriptor == AKAZE::DESCRIPTOR_MLDB_UPRIGHT) { + int num_keypoints = (int)kpts.size(); + + // Group keypoints by their level to process each level independently + std::vector> kpts_by_level(uPyr.size()); + for (int i = 0; i < num_keypoints; i++) { + int level = kpts[i].class_id; + if (level >= 0 && level < (int)uPyr.size()) { + kpts_by_level[level].push_back(i); + } + } + + // Process each level independently on GPU + for (size_t level = 0; level < uPyr.size(); level++) { + if (kpts_by_level[level].empty()) + continue; + + const Evolution& e = uPyr[level]; + int num_kpts_level = (int)kpts_by_level[level].size(); + + // Create GPU buffers for keypoints in this level + std::vector kpts_x, kpts_y, kpts_size, kpts_angle; + + for (int idx : kpts_by_level[level]) { + kpts_x.push_back(kpts[idx].pt.x); + kpts_y.push_back(kpts[idx].pt.y); + kpts_size.push_back(kpts[idx].size); + kpts_angle.push_back(kpts[idx].angle); + } + + Mat keypoints_x_mat(num_kpts_level, 1, CV_32F, kpts_x.data()); + Mat keypoints_y_mat(num_kpts_level, 1, CV_32F, kpts_y.data()); + Mat keypoints_size_mat(num_kpts_level, 1, CV_32F, kpts_size.data()); + Mat keypoints_angle_mat(num_kpts_level, 1, CV_32F, kpts_angle.data()); + + UMat keypoints_x, keypoints_y, keypoints_size, keypoints_angle; + keypoints_x_mat.copyTo(keypoints_x); + keypoints_y_mat.copyTo(keypoints_y); + keypoints_size_mat.copyTo(keypoints_size); + keypoints_angle_mat.copyTo(keypoints_angle); + + // Get level data directly without concatenation + UMat Lx_level = e.Lx; + UMat Ly_level = e.Ly; + UMat Lt_level = e.Lt; + + // Allocate descriptor buffer for this level + UMat descriptors_level(num_kpts_level, descriptor_size, CV_8UC1); + + // Compute descriptors for this level on GPU + ocl_compute_mldb_descriptor_level(keypoints_x, keypoints_y, keypoints_size, keypoints_angle, + Lx_level, Ly_level, Lt_level, + e.Lx.rows, e.Lx.cols, e.octave_ratio, + descriptor_size, descriptors_level); + + // Copy descriptors to final output + Mat descriptors_level_mat; + descriptors_level.copyTo(descriptors_level_mat); + + Mat descriptors_mat = desc.getMat(); + for (int i = 0; i < num_kpts_level; i++) { + int global_idx = kpts_by_level[level][i]; + for (int j = 0; j < descriptor_size; j++) { + descriptors_mat.at(global_idx, j) = descriptors_level_mat.at(i, j); + } + } + } + } else { + // For non-MLDB descriptors, fall back to CPU implementation + Mat desc_mat = desc.getMat(); + Compute_Descriptors(kpts, desc_mat); + } +} +#endif // HAVE_OPENCL + /** * @brief This method finds extrema in the nonlinear scale space * @param keypoints_by_layers Output vectors of detected keypoints; one vector for each evolution level diff --git a/modules/features2d/src/kaze/AKAZEFeatures.h b/modules/features2d/src/kaze/AKAZEFeatures.h index 512f553886..7e78f5b9c7 100644 --- a/modules/features2d/src/kaze/AKAZEFeatures.h +++ b/modules/features2d/src/kaze/AKAZEFeatures.h @@ -94,6 +94,12 @@ private: void Do_Subpixel_Refinement(std::vector& keypoints_by_layers, std::vector& kpts); +public: + void GetEvolutionPyramid(UMatPyramid& uPyr); + void Create_Nonlinear_Scale_Space_UMat(InputArray img, UMatPyramid& uPyr); + void Feature_Detection_UMat(UMatPyramid& uPyr, std::vector& kpts); + void Compute_Descriptors_UMat(std::vector& kpts, OutputArray desc, UMatPyramid& uPyr); + /// Feature description methods void Compute_Keypoints_Orientation(std::vector& kpts) const; diff --git a/modules/features2d/src/kaze/KAZEFeatures.cpp b/modules/features2d/src/kaze/KAZEFeatures.cpp index 30e0390359..8d27df041a 100644 --- a/modules/features2d/src/kaze/KAZEFeatures.cpp +++ b/modules/features2d/src/kaze/KAZEFeatures.cpp @@ -23,6 +23,8 @@ #include "../precomp.hpp" #include "KAZEFeatures.h" #include "utils.h" +#include "opencl_kernels_features2d.hpp" +#include "akaze_diffusion.hpp" namespace cv { @@ -196,13 +198,13 @@ void KAZEFeatures::Feature_Detection(std::vector& kpts) class MultiscaleDerivativesKAZEInvoker : public ParallelLoopBody { public: - explicit MultiscaleDerivativesKAZEInvoker(std::vector& ev) : evolution_(&ev) + explicit MultiscaleDerivativesKAZEInvoker(TPyramid& ev) : evolution_(&ev) { } void operator()(const Range& range) const CV_OVERRIDE { - std::vector& evolution = *evolution_; + TPyramid& evolution = *evolution_; for (int i = range.start; i < range.end; i++) { compute_scharr_derivatives(evolution[i].Lsmooth, evolution[i].Lx, 1, 0, evolution[i].sigma_size); @@ -238,14 +240,14 @@ void KAZEFeatures::Compute_Multiscale_Derivatives(void) class FindExtremumKAZEInvoker : public ParallelLoopBody { public: - explicit FindExtremumKAZEInvoker(std::vector& ev, std::vector >& kpts_par, + explicit FindExtremumKAZEInvoker(TPyramid& ev, std::vector >& kpts_par, const KAZEOptions& options) : evolution_(&ev), kpts_par_(&kpts_par), options_(options) { } void operator()(const Range& range) const CV_OVERRIDE { - std::vector& evolution = *evolution_; + TPyramid& evolution = *evolution_; std::vector >& kpts_par = *kpts_par_; for (int i = range.start; i < range.end; i++) { @@ -300,7 +302,7 @@ public: } private: - std::vector* evolution_; + TPyramid* evolution_; std::vector >* kpts_par_; KAZEOptions options_; }; @@ -499,7 +501,7 @@ void KAZEFeatures::Do_Subpixel_Refinement(std::vector &kpts) { class KAZE_Descriptor_Invoker : public ParallelLoopBody { public: - KAZE_Descriptor_Invoker(std::vector &kpts, Mat &desc, std::vector& evolution, const KAZEOptions& options) + KAZE_Descriptor_Invoker(std::vector &kpts, Mat &desc, TPyramid& evolution, const KAZEOptions& options) : kpts_(&kpts) , desc_(&desc) , evolution_(&evolution) @@ -515,7 +517,7 @@ public: { std::vector &kpts = *kpts_; Mat &desc = *desc_; - std::vector &evolution = *evolution_; + TPyramid &evolution = *evolution_; for (int i = range.start; i < range.end; i++) { @@ -547,7 +549,7 @@ private: std::vector * kpts_; Mat * desc_; - std::vector * evolution_; + TPyramid * evolution_; KAZEOptions options_; }; @@ -582,7 +584,7 @@ void KAZEFeatures::Feature_Description(std::vector &kpts, Mat &desc) * @note The orientation is computed using a similar approach as described in the * original SURF method. See Bay et al., Speeded Up Robust Features, ECCV 2006 */ -void KAZEFeatures::Compute_Main_Orientation(KeyPoint &kpt, const std::vector& evolution_, const KAZEOptions& options) +void KAZEFeatures::Compute_Main_Orientation(KeyPoint &kpt, const TPyramid& evolution_, const KAZEOptions& options) { int ix = 0, iy = 0, idx = 0, s = 0, level = 0; float xf = 0.0, yf = 0.0, gweight = 0.0; @@ -671,7 +673,7 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Upright_Descriptor_64(const KeyPoint &kpt float fx = 0.0, fy = 0.0, res1 = 0.0, res2 = 0.0, res3 = 0.0, res4 = 0.0; int dsize = 0, scale = 0, level = 0; - std::vector& evolution = *evolution_; + TPyramid& evolution = *evolution_; // Subregion centers for the 4x4 gaussian weighting float cx = -0.5f, cy = 0.5f; @@ -799,7 +801,7 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Descriptor_64(const KeyPoint &kpt, float int kx = 0, ky = 0, i = 0, j = 0, dcount = 0; int dsize = 0, scale = 0, level = 0; - std::vector& evolution = *evolution_; + TPyramid& evolution = *evolution_; // Subregion centers for the 4x4 gaussian weighting float cx = -0.5f, cy = 0.5f; @@ -933,7 +935,7 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Upright_Descriptor_128(const KeyPoint &kp // Subregion centers for the 4x4 gaussian weighting float cx = -0.5f, cy = 0.5f; - std::vector& evolution = *evolution_; + TPyramid& evolution = *evolution_; // Set the descriptor size and the sample and pattern sizes dsize = 128; @@ -1082,7 +1084,7 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Descriptor_128(const KeyPoint &kpt, float int kx = 0, ky = 0, i = 0, j = 0, dcount = 0; int dsize = 0, scale = 0, level = 0; - std::vector& evolution = *evolution_; + TPyramid& evolution = *evolution_; // Subregion centers for the 4x4 gaussian weighting float cx = -0.5f, cy = 0.5f; @@ -1217,4 +1219,254 @@ void KAZE_Descriptor_Invoker::Get_KAZE_Descriptor_128(const KeyPoint &kpt, float } } +#ifdef HAVE_OPENCL + +// Compute GaussianBlur kernel size matching KAZE's gaussian_2D_convolution formula +static inline int kaze_ksize_for_sigma(float sigma) +{ + return cvCeil(2.0f * (1.0f + (sigma - 0.8f) / 0.3f)) | 1; +} + +// Dispatch the KAZE upright descriptor kernel for one pyramid level. +// lx_step: row stride of Lx/Ly in float elements (= Lx.cols for continuous mat). +static inline bool +ocl_compute_kaze_upright_descriptor_level(InputArray Lx_, InputArray Ly_, + InputArray keypoints_x_, InputArray keypoints_y_, + InputArray keypoints_size_, + OutputArray descriptors_, bool extended) +{ + UMat Lx = Lx_.getUMat(); + UMat Ly = Ly_.getUMat(); + UMat keypoints_x = keypoints_x_.getUMat(); + UMat keypoints_y = keypoints_y_.getUMat(); + UMat keypoints_size = keypoints_size_.getUMat(); + UMat descriptors = descriptors_.getUMat(); + + int nkeypoints = (int)keypoints_x.total(); + // For continuous UMat, step_in_float_elements == cols + int lx_step = Lx.cols; + int lx_rows = Lx.rows; + int lx_cols = Lx.cols; + + size_t globalSize[] = {(size_t)nkeypoints}; + + const char* kernel_name = extended ? "KAZE_compute_upright_descriptor_128" + : "KAZE_compute_upright_descriptor_64"; + ocl::Kernel ker(kernel_name, ocl::features2d::kaze_oclsrc); + if (ker.empty()) + return false; + + return ker.args( + ocl::KernelArg::PtrReadOnly(Lx), + ocl::KernelArg::PtrReadOnly(Ly), + lx_step, lx_rows, lx_cols, + ocl::KernelArg::PtrReadOnly(keypoints_x), + ocl::KernelArg::PtrReadOnly(keypoints_y), + ocl::KernelArg::PtrReadOnly(keypoints_size), + ocl::KernelArg::PtrWriteOnly(descriptors), + nkeypoints).run(1, globalSize, 0, false); +} + +void KAZEFeatures::Create_Nonlinear_Scale_Space_UMat(InputArray image, UTPyramid& uPyr) +{ + CV_Assert(evolution_.size() > 0); + + // Initialize UTPyramid metadata from the pre-built CPU evolution + uPyr.resize(evolution_.size()); + for (size_t i = 0; i < evolution_.size(); i++) + { + uPyr[i].etime = evolution_[i].etime; + uPyr[i].esigma = evolution_[i].esigma; + uPyr[i].octave = evolution_[i].octave; + uPyr[i].sublevel = evolution_[i].sublevel; + uPyr[i].sigma_size = evolution_[i].sigma_size; + } + + // Level 0: soffset blur (-> Lt), then sderivatives blur (-> Lsmooth) + { + int k0 = kaze_ksize_for_sigma(options_.soffset); + int k1 = kaze_ksize_for_sigma(options_.sderivatives); + GaussianBlur(image, uPyr[0].Lt, Size(k0, k0), options_.soffset, options_.soffset, BORDER_REPLICATE); + GaussianBlur(uPyr[0].Lt, uPyr[0].Lsmooth, Size(k1, k1), options_.sderivatives, options_.sderivatives, BORDER_REPLICATE); + } + + // Compute kcontrast on CPU using first-level Lt (one GPU→CPU sync) + { + Mat lt_cpu; + uPyr[0].Lt.copyTo(lt_cpu); + options_.kcontrast = compute_k_percentile(lt_cpu, options_.kcontrast_percentille, + options_.sderivatives, options_.kcontrast_bins, 0, 0); + } + + // Build scale space (levels 1 .. n-1) + { + int k1 = kaze_ksize_for_sigma(options_.sderivatives); + UMat Lx_diff, Ly_diff, Lflow, Lstep; + + for (size_t i = 1; i < uPyr.size(); i++) + { + // Copy diffused previous level and compute Lsmooth from it + uPyr[i - 1].Lt.copyTo(uPyr[i].Lt); + GaussianBlur(uPyr[i - 1].Lt, uPyr[i].Lsmooth, + Size(k1, k1), options_.sderivatives, options_.sderivatives, BORDER_REPLICATE); + + // Scale=1 Scharr for diffusivity (not stored; used only to drive diffusion) + Scharr(uPyr[i].Lsmooth, Lx_diff, CV_32F, 1, 0, 1.0, 0.0, BORDER_DEFAULT); + Scharr(uPyr[i].Lsmooth, Ly_diff, CV_32F, 0, 1, 1.0, 0.0, BORDER_DEFAULT); + + compute_diffusivity(Lx_diff, Ly_diff, Lflow, options_.kcontrast, options_.diffusivity); + + // Fast Explicit Diffusion steps + const std::vector& tsteps = tsteps_[i - 1]; + for (size_t j = 0; j < tsteps.size(); j++) + { + non_linear_diffusion_step(uPyr[i].Lt, Lflow, Lstep, tsteps[j] * 0.5f); + add(uPyr[i].Lt, Lstep, uPyr[i].Lt); + } + } + } + + // Compute sigma_size-scaled multiscale derivatives and Ldet for all levels. + // Matches CPU Compute_Multiscale_Derivatives + Compute_Detector_Response exactly: + // Lx = s * scharr_s(Lsmooth) + // Ly = s * scharr_s(Lsmooth) + // Lxx = s^2 * scharr_s(Lx_unscaled) [applied before scaling Lx] + // Lyy = s^2 * scharr_s(Ly_unscaled) + // Lxy = s^2 * scharr_s_y(Lx_unscaled) + // Ldet = Lxx*Lyy - Lxy^2 + for (size_t i = 0; i < uPyr.size(); i++) + { + int s = uPyr[i].sigma_size; + + Mat kx_x, ky_x, kx_y, ky_y; + compute_derivative_kernels(kx_x, ky_x, 1, 0, s); + compute_derivative_kernels(kx_y, ky_y, 0, 1, s); + + // First-order (unscaled) from Lsmooth + UMat Lx_unscaled, Ly_unscaled; + sepFilter2D(uPyr[i].Lsmooth, Lx_unscaled, CV_32F, kx_x, ky_x); + sepFilter2D(uPyr[i].Lsmooth, Ly_unscaled, CV_32F, kx_y, ky_y); + + // Second-order from unscaled first-order (must happen before scaling Lx/Ly) + UMat Lxx_unscaled, Lxy_unscaled, Lyy_unscaled; + sepFilter2D(Lx_unscaled, Lxx_unscaled, CV_32F, kx_x, ky_x); + sepFilter2D(Ly_unscaled, Lyy_unscaled, CV_32F, kx_y, ky_y); + sepFilter2D(Lx_unscaled, Lxy_unscaled, CV_32F, kx_y, ky_y); + + // Scale first-order and store in pyramid (used by descriptor kernel) + Lx_unscaled.convertTo(uPyr[i].Lx, CV_32F, (double)s); + Ly_unscaled.convertTo(uPyr[i].Ly, CV_32F, (double)s); + + // Scale second-order and compute Ldet + double s2 = (double)s * s; + UMat Lxx, Lxy, Lyy; + Lxx_unscaled.convertTo(Lxx, CV_32F, s2); + Lxy_unscaled.convertTo(Lxy, CV_32F, s2); + Lyy_unscaled.convertTo(Lyy, CV_32F, s2); + + UMat tmp1, tmp2; + multiply(Lxx, Lyy, tmp1); + multiply(Lxy, Lxy, tmp2); + subtract(tmp1, tmp2, uPyr[i].Ldet); + } +} + +void KAZEFeatures::Feature_Detection_UMat(UTPyramid& uPyr, std::vector& kpts) +{ + // GPU built the scale space; download Ldet for all levels, then run CPU detection. + // The CPU Determinant_Hessian + Do_Subpixel_Refinement only need Ldet. + for (size_t i = 0; i < uPyr.size(); i++) + uPyr[i].Ldet.copyTo(evolution_[i].Ldet); // single batch sync on first copyTo + + Determinant_Hessian(kpts); + Do_Subpixel_Refinement(kpts); +} + +void KAZEFeatures::Compute_Descriptors_UMat(std::vector& kpts, OutputArray desc, + UTPyramid& uPyr) +{ + if (kpts.empty()) + { + desc.release(); + return; + } + + int descriptor_size = options_.extended ? 128 : 64; + + // Group keypoint indices by pyramid level + std::vector> kpts_by_level(uPyr.size()); + for (size_t i = 0; i < kpts.size(); i++) + { + int lv = kpts[i].class_id; + if (lv >= 0 && lv < (int)uPyr.size()) + kpts_by_level[lv].push_back((int)i); + } + + // Pre-allocate output matrix; we'll write into row-ranges per level + desc.create((int)kpts.size(), descriptor_size, CV_32F); + Mat out_mat = desc.getMat(); + + // Track the keypoint ordering as we iterate levels + // (kpts_by_level iterates in ascending level order, matching kpts order after sorting) + // We use a small vector to map global keypoint index -> descriptor row + std::vector kpt_to_row(kpts.size(), -1); + { + // Determine final output row for each keypoint based on iteration order + int row = 0; + for (size_t lv = 0; lv < uPyr.size(); lv++) + for (int ki : kpts_by_level[lv]) + kpt_to_row[ki] = row++; + } + + for (size_t lv = 0; lv < uPyr.size(); lv++) + { + const std::vector& indices = kpts_by_level[lv]; + int n = (int)indices.size(); + if (n == 0) + continue; + + // Pack keypoint data into UMats for this level + UMat kpts_x(n, 1, CV_32F); + UMat kpts_y(n, 1, CV_32F); + UMat kpts_sz(n, 1, CV_32F); + { + Mat mx = kpts_x.getMat(ACCESS_WRITE); + Mat my = kpts_y.getMat(ACCESS_WRITE); + Mat ms = kpts_sz.getMat(ACCESS_WRITE); + for (int i = 0; i < n; i++) + { + const KeyPoint& kp = kpts[indices[i]]; + mx.at(i) = kp.pt.x; + my.at(i) = kp.pt.y; + ms.at(i) = kp.size; + } + } + + // Compute descriptors for this level on GPU + UMat desc_level(n, descriptor_size, CV_32F); + bool ok = ocl_compute_kaze_upright_descriptor_level( + uPyr[lv].Lx, uPyr[lv].Ly, + kpts_x, kpts_y, kpts_sz, + desc_level, options_.extended); + + if (!ok) + { + // GPU failed — fall back to full CPU path + Mat cpu_desc; + Feature_Description(kpts, cpu_desc); + cpu_desc.copyTo(desc); + return; + } + + // Download and copy into correct rows of the output matrix + Mat level_mat = desc_level.getMat(ACCESS_READ); + for (int i = 0; i < n; i++) + { + int row = kpt_to_row[indices[i]]; + level_mat.row(i).copyTo(out_mat.row(row)); + } + } +} +#endif // HAVE_OPENCL + } diff --git a/modules/features2d/src/kaze/KAZEFeatures.h b/modules/features2d/src/kaze/KAZEFeatures.h index 934dfff29d..3fdd4d7cd5 100644 --- a/modules/features2d/src/kaze/KAZEFeatures.h +++ b/modules/features2d/src/kaze/KAZEFeatures.h @@ -28,7 +28,7 @@ private: /// Parameters of the Nonlinear diffusion class KAZEOptions options_; ///< Configuration options for KAZE - std::vector evolution_; ///< Vector of nonlinear diffusion evolution + TPyramid evolution_; ///< Vector of nonlinear diffusion evolution /// Vector of keypoint vectors for finding extrema in multiple threads std::vector > kpts_par_; @@ -49,7 +49,13 @@ public: int Create_Nonlinear_Scale_Space(const cv::Mat& img); void Feature_Detection(std::vector& kpts); void Feature_Description(std::vector& kpts, cv::Mat& desc); - static void Compute_Main_Orientation(cv::KeyPoint& kpt, const std::vector& evolution_, const KAZEOptions& options); + static void Compute_Main_Orientation(cv::KeyPoint& kpt, const TPyramid& evolution_, const KAZEOptions& options); + +#ifdef HAVE_OPENCL + void Create_Nonlinear_Scale_Space_UMat(InputArray img, UTPyramid& uPyr); + void Feature_Detection_UMat(UTPyramid& uPyr, std::vector& kpts); + void Compute_Descriptors_UMat(std::vector& kpts, OutputArray desc, UTPyramid& uPyr); +#endif /// Feature Detection Methods void Compute_KContrast(const cv::Mat& img, const float& kper); diff --git a/modules/features2d/src/kaze/TEvolution.h b/modules/features2d/src/kaze/TEvolution.h index 36052b3e64..df08afd4d3 100644 --- a/modules/features2d/src/kaze/TEvolution.h +++ b/modules/features2d/src/kaze/TEvolution.h @@ -13,9 +13,10 @@ namespace cv /* ************************************************************************* */ /// KAZE/A-KAZE nonlinear diffusion filtering evolution -struct TEvolution +template +struct TEvolutionT { - TEvolution() { + TEvolutionT() { etime = 0.0f; esigma = 0.0f; octave = 0; @@ -23,11 +24,11 @@ struct TEvolution sigma_size = 0; } - Mat Lx, Ly; ///< First order spatial derivatives - Mat Lxx, Lxy, Lyy; ///< Second order spatial derivatives - Mat Lt; ///< Evolution image - Mat Lsmooth; ///< Smoothed image - Mat Ldet; ///< Detector response + MatType Lx, Ly; ///< First order spatial derivatives + MatType Lxx, Lxy, Lyy; ///< Second order spatial derivatives + MatType Lt; ///< Evolution image + MatType Lsmooth; ///< Smoothed image + MatType Ldet; ///< Detector response float etime; ///< Evolution time float esigma; ///< Evolution sigma. For linear diffusion t = sigma^2 / 2 @@ -36,6 +37,11 @@ struct TEvolution int sigma_size; ///< Integer esigma. For computing the feature detector responses }; +typedef TEvolutionT TEvolution; +typedef TEvolutionT UTEvolution; +typedef std::vector TPyramid; +typedef std::vector UTPyramid; + } #endif diff --git a/modules/features2d/src/kaze/akaze_diffusion.hpp b/modules/features2d/src/kaze/akaze_diffusion.hpp new file mode 100644 index 0000000000..d274c1bb7d --- /dev/null +++ b/modules/features2d/src/kaze/akaze_diffusion.hpp @@ -0,0 +1,218 @@ +// Internal header: shared nonlinear diffusion utilities used by both +// AKAZEFeatures.cpp and KAZEFeatures.cpp. +// Not part of the public API. + +#ifndef __OPENCV_FEATURES_2D_AKAZE_DIFFUSION_HPP__ +#define __OPENCV_FEATURES_2D_AKAZE_DIFFUSION_HPP__ + +#include "../precomp.hpp" +#include "nldiffusion_functions.h" +#ifdef HAVE_OPENCL +#include "opencl_kernels_features2d.hpp" +#endif + +namespace cv +{ + +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(); + + 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; + + if (row == 0) { + lt_c = Lt.ptr(0) + 1; + lf_c = Lf.ptr(0) + 1; + lt_b = Lt.ptr(1) + 1; + lf_b = Lf.ptr(1) + 1; + + dst = Lstep.ptr(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; + } + + dst[cols] = 0.0f; + ++row; + } + + int middle_end = std::min(Lt.rows - 1, row_end); + for (; row < middle_end; ++row) + { + lt_a = Lt.ptr(row - 1); + lf_a = Lf.ptr(row - 1); + lt_c = Lt.ptr(row ); + lf_c = Lf.ptr(row ); + lt_b = Lt.ptr(row + 1); + lf_b = Lf.ptr(row + 1); + dst = Lstep.ptr(row); + + 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++; + + 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; + } + + 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; + } + + if (row_end == Lt.rows) { + lt_a = Lt.ptr(row - 1) + 1; + lf_a = Lf.ptr(row - 1) + 1; + lt_c = Lt.ptr(row ) + 1; + lf_c = Lf.ptr(row ) + 1; + + dst = Lstep.ptr(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; + } + + 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 CV_OVERRIDE + { + 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(InputArray Lt_, InputArray Lf_, OutputArray Lstep_, float step_size) +{ + if(!Lt_.isContinuous()) + return false; + + UMat Lt = Lt_.getUMat(); + UMat Lf = Lf_.getUMat(); + UMat Lstep = Lstep_.getUMat(); + + 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, false); +} + +static inline bool +ocl_pm_g2(InputArray Lx_, InputArray Ly_, OutputArray Lflow_, float kcontrast) +{ + UMat Lx = Lx_.getUMat(); + UMat Ly = Ly_.getUMat(); + UMat Lflow = Lflow_.getUMat(); + + 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, false); +} +#endif // HAVE_OPENCL + +static inline void +non_linear_diffusion_step(InputArray Lt_, InputArray Lf_, OutputArray Lstep_, float step_size) +{ + CV_INSTRUMENT_REGION(); + + Lstep_.create(Lt_.size(), Lt_.type()); + + CV_OCL_RUN(Lt_.isUMat() && Lf_.isUMat() && Lstep_.isUMat(), + ocl_non_linear_diffusion_step(Lt_, Lf_, Lstep_, step_size)); + + Mat Lt = Lt_.getMat(); + Mat Lf = Lf_.getMat(); + Mat Lstep = Lstep_.getMat(); + parallel_for_(Range(0, Lt.rows), NonLinearScalarDiffusionStep(Lt, Lf, Lstep, step_size)); +} + +static inline void +compute_diffusivity(InputArray Lx, InputArray Ly, OutputArray Lflow, float kcontrast, KAZE::DiffusivityType 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(Lx.isUMat() && Ly.isUMat() && Lflow.isUMat(), 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_(Error::StsError, ("Diffusivity is not supported: %d", static_cast(diffusivity))); + break; + } +} + +} // namespace cv + +#endif // __OPENCV_FEATURES_2D_AKAZE_DIFFUSION_HPP__ diff --git a/modules/features2d/src/opencl/akaze.cl b/modules/features2d/src/opencl/akaze.cl index 40f5c2b59c..7664a8752b 100644 --- a/modules/features2d/src/opencl/akaze.cl +++ b/modules/features2d/src/opencl/akaze.cl @@ -120,3 +120,527 @@ AKAZE_compute_determinant(__global const float* lxx, __global const float* lxy, dst[i] = (lxx[i] * lyy[i] - lxy[i] * lxy[i]) * sigma; } + +/** + * @brief Find scale space extrema in 3x3 neighborhood + * @details Detects local maxima in Hessian response within 3x3 neighborhood + * + * @param ldet Hessian determinant response + * @param rows image height + * @param cols image width + * @param threshold detector response threshold + * @param border border to ignore + * @param keypoint_mask output binary mask of detected keypoints + */ +__kernel void +AKAZE_find_extrema_same_scale(__global const float* ldet, int rows, int cols, + float threshold, int border, __global uchar* keypoint_mask) +{ + int x = get_global_id(0); + int y = get_global_id(1); + + // Check bounds + if (x < border || x >= cols - border || y < border || y >= rows - border) + { + return; + } + + int idx = y * cols + x; + float value = ldet[idx]; + + // Filter by threshold + if (value <= threshold) + { + keypoint_mask[idx] = 0; + return; + } + + // 3x3 local maxima check + // Check horizontal neighbors + if (value <= ldet[idx - 1] || value <= ldet[idx + 1]) + { + keypoint_mask[idx] = 0; + return; + } + + // Check vertical neighbors + int prev_row = (y - 1) * cols + x; + int next_row = (y + 1) * cols + x; + if (value <= ldet[prev_row] || value <= ldet[next_row]) + { + keypoint_mask[idx] = 0; + return; + } + + // Check diagonal neighbors + if (value <= ldet[prev_row - 1] || value <= ldet[prev_row + 1] || + value <= ldet[next_row - 1] || value <= ldet[next_row + 1]) + { + keypoint_mask[idx] = 0; + return; + } + + // This is a local maximum + keypoint_mask[idx] = 1; +} + +/** + * @brief Cross-scale non-maximum suppression (lower scale filtering) + * @details For each keypoint in current level, project to lower level and suppress if weaker + * + * @param keypoints_current keypoint mask for current scale level + * @param keypoints_lower keypoint mask for lower scale level + * @param ldet_current Hessian response for current level + * @param ldet_lower Hessian response for lower level + * @param rows image height + * @param cols image width + * @param diff_ratio ratio to project from current to lower level + * @param search_radius search radius in lower level + */ +__kernel void +AKAZE_cross_scale_filter_lower(__global const uchar* keypoints_current, + __global uchar* keypoints_lower, + __global const float* ldet_current, + __global const float* ldet_lower, + int rows, int cols, + int lower_rows, int lower_cols, + int diff_ratio, int search_radius) +{ + int x = get_global_id(0); + int y = get_global_id(1); + + if (x >= cols || y >= rows) + return; + + int idx = y * cols + x; + + // Only process keypoints in current level + if (keypoints_current[idx] == 0) + return; + + // Project to lower scale level + int p_x = x * diff_ratio; + int p_y = y * diff_ratio; + + // Search for neighbor in lower level within radius + int radius_sq = search_radius * search_radius; + int found = 0; + int neighbor_idx = -1; + + // Brute force search within radius + for (int dy = -search_radius; dy <= search_radius; dy++) + { + for (int dx = -search_radius; dx <= search_radius; dx++) + { + int nx = p_x + dx; + int ny = p_y + dy; + + if (nx >= 0 && nx < lower_cols && ny >= 0 && ny < lower_rows) + { + int nidx = ny * lower_cols + nx; + if (keypoints_lower[nidx] == 1) + { + if (dx * dx + dy * dy <= radius_sq) + { + found = 1; + neighbor_idx = nidx; + break; + } + } + } + } + if (found) break; + } + + // If neighbor found and current response is higher, suppress neighbor + if (found && ldet_current[idx] > ldet_lower[neighbor_idx]) + { + keypoints_lower[neighbor_idx] = 0; + } +} + +/** + * @brief Cross-scale non-maximum suppression (upper scale filtering) + * @details For each keypoint in current level, project to upper level and suppress if weaker + * + * @param keypoints_current keypoint mask for current scale level + * @param keypoints_upper keypoint mask for upper scale level + * @param ldet_current Hessian response for current level + * @param ldet_upper Hessian response for upper level + * @param rows image height + * @param cols image width + * @param diff_ratio ratio to project from current to upper level + * @param search_radius search radius in upper level + */ +__kernel void +AKAZE_cross_scale_filter_upper(__global const uchar* keypoints_current, + __global uchar* keypoints_upper, + __global const float* ldet_current, + __global const float* ldet_upper, + int rows, int cols, + int upper_rows, int upper_cols, + int diff_ratio, int search_radius) +{ + int x = get_global_id(0); + int y = get_global_id(1); + + if (x >= cols || y >= rows) + return; + + int idx = y * cols + x; + + // Only process keypoints in current level + if (keypoints_current[idx] == 0) + return; + + // Project to upper scale level + int p_x = x / diff_ratio; + int p_y = y / diff_ratio; + + // Search for neighbor in upper level within radius + int radius_sq = search_radius * search_radius; + int found = 0; + int neighbor_idx = -1; + + // Brute force search within radius + for (int dy = -search_radius; dy <= search_radius; dy++) + { + for (int dx = -search_radius; dx <= search_radius; dx++) + { + int nx = p_x + dx; + int ny = p_y + dy; + + if (nx >= 0 && nx < upper_cols && ny >= 0 && ny < upper_rows) + { + int nidx = ny * upper_cols + nx; + if (keypoints_upper[nidx] == 1) + { + if (dx * dx + dy * dy <= radius_sq) + { + found = 1; + neighbor_idx = nidx; + break; + } + } + } + } + if (found) break; + } + + // If neighbor found and current response is higher, suppress neighbor + if (found && ldet_current[idx] > ldet_upper[neighbor_idx]) + { + keypoints_upper[neighbor_idx] = 0; + } +} + +/** + * @brief Combined subpixel refinement and orientation for single level + * @details Processes subpixel refinement and orientation together to eliminate intermediate transfers + * + * @param keypoints keypoint mask for this level + * @param ldet Hessian response for this level + * @param Lx gradient in x direction for this level + * @param Ly gradient in y direction for this level + * @param rows image height + * @param cols image width + * @param octave_ratio scale ratio for this level + * @param esigma evolution sigma for this level + * @param octave octave number for this level + * @param level evolution level index + * @param output_count atomic counter for number of refined keypoints + * @param output_x output x coordinates + * @param output_y output y coordinates + * @param output_size output keypoint sizes + * @param output_response output keypoint responses + * @param output_octave output octave numbers + * @param output_class_id output level indices + * @param output_angle output orientation angles + * @param max_output maximum number of keypoints that can be output + */ +__kernel void +AKAZE_subpixel_refinement_orientation(__global const uchar* keypoints, + __global const float* ldet, + __global const float* Lx, + __global const float* Ly, + int rows, int cols, + float octave_ratio, float esigma, int octave, int level, + __global int* output_count, + __global float* output_x, + __global float* output_y, + __global float* output_size, + __global float* output_response, + __global int* output_octave, + __global int* output_class_id, + __global float* output_angle, + int max_output) +{ + int x = get_global_id(0); + int y = get_global_id(1); + + if (x >= cols || y >= rows) + return; + + int idx = y * cols + x; + + // Only process keypoints + if (keypoints[idx] == 0) + return; + + // Compute gradient for subpixel refinement + float Dx = 0.5f * (ldet[idx + 1] - ldet[idx - 1]); + float Dy = 0.5f * (ldet[idx + cols] - ldet[idx - cols]); + + // Compute Hessian + float Dxx = ldet[idx + 1] + ldet[idx - 1] - 2.0f * ldet[idx]; + float Dyy = ldet[idx + cols] + ldet[idx - cols] - 2.0f * ldet[idx]; + float Dxy = 0.25f * (ldet[idx + cols + 1] + ldet[idx - cols - 1] - + ldet[idx - cols + 1] - ldet[idx + cols - 1]); + + // Solve 2x2 linear system + float det = Dxx * Dyy - Dxy * Dxy; + + if (fabs(det) < 1e-10f) + return; + + float inv_det = 1.0f / det; + float dx = (-Dyy * Dx + Dxy * Dy) * inv_det; + float dy = (Dxy * Dx - Dxx * Dy) * inv_det; + + if (fabs(dx) > 1.0f || fabs(dy) > 1.0f) + return; + + int out_idx = atomic_inc(output_count); + + if (out_idx >= max_output) + return; + + float refined_x = x * octave_ratio + dx * octave_ratio + 0.5f * (octave_ratio - 1.0f); + float refined_y = y * octave_ratio + dy * octave_ratio + 0.5f * (octave_ratio - 1.0f); + + output_x[out_idx] = refined_x; + output_y[out_idx] = refined_y; + output_size[out_idx] = esigma * 3.0f; // derivative_factor(1.5) * 2.0(diameter) = 3.0 + output_response[out_idx] = ldet[idx]; + output_octave[out_idx] = octave; + output_class_id[out_idx] = level; + + // Compute orientation using gradient histogram + float scale_f = 0.5f * output_size[out_idx] / octave_ratio; + int scale = (int)(scale_f + 0.5f); + float x0_f = refined_x / octave_ratio; + int x0 = (int)(x0_f + (x0_f >= 0 ? 0.5f : -0.5f)); + float y0_f = refined_y / octave_ratio; + int y0 = (int)(y0_f + (y0_f >= 0 ? 0.5f : -0.5f)); + + float histogram[36]; + for (int b = 0; b < 36; b++) + histogram[b] = 0.0f; + + const int radius = 6 * scale; + const int radius_sq = radius * radius; + + for (int dy = -radius; dy <= radius; dy++) { + for (int dx = -radius; dx <= radius; dx++) { + if (dx * dx + dy * dy > radius_sq) + continue; + + int sx = x0 + dx; + int sy = y0 + dy; + + if (sx < 0 || sx >= cols || sy < 0 || sy >= rows) + continue; + + int sidx = sy * cols + sx; + float gx = Lx[sidx]; + float gy = Ly[sidx]; + + float magnitude = sqrt(gx * gx + gy * gy); + float angle = atan2(gy, gx); + if (angle < 0.0f) + angle += 2.0f * M_PI_F; + + float dist = sqrt((float)(dx * dx + dy * dy)); + float weight = exp(-dist * dist / (2.0f * scale * scale)); + + int bin = (int)(angle * 36 / (2.0f * M_PI_F)); + bin = bin % 36; + histogram[bin] += magnitude * weight; + } + } + + float max_hist = 0.0f; + int max_bin = 0; + for (int b = 0; b < 36; b++) { + if (histogram[b] > max_hist) { + max_hist = histogram[b]; + max_bin = b; + } + } + + if (max_hist == 0.0f) + { + output_angle[out_idx] = 0.0f; + return; + } + + int prev_bin = (max_bin - 1 + 36) % 36; + int next_bin = (max_bin + 1) % 36; + + float prev_val = histogram[prev_bin]; + float curr_val = histogram[max_bin]; + float next_val = histogram[next_bin]; + + float denom = prev_val - 2.0f * curr_val + next_val; + float delta = (fabs(denom) > 1e-10f) ? 0.5f * (prev_val - next_val) / denom : 0.0f; + float refined_angle = (max_bin + delta) * (2.0f * M_PI_F / 36); + + output_angle[out_idx] = refined_angle * 180.0f / M_PI_F; +} + +/** + * @brief Compute MLDB descriptor for a keypoint (single level, upright version) + * @details Samples binary comparisons from gradient responses matching CPU implementation + * + * @param keypoints_x keypoint x coordinates + * @param keypoints_y keypoint y coordinates + * @param keypoints_size keypoint sizes + * @param keypoints_angle keypoint orientations (unused in upright version) + * @param num_keypoints number of keypoints + * @param Lx gradient in x direction + * @param Ly gradient in y direction + * @param Lt flow response + * @param rows image height + * @param cols image width + * @param octave_ratio octave ratio for this level + * @param descriptor_size size of descriptor in bytes + * @param output_descriptors output descriptor matrix + */ +__kernel void +AKAZE_compute_mldb_descriptor_level(__global const float* keypoints_x, + __global const float* keypoints_y, + __global const float* keypoints_size, + __global const float* keypoints_angle, + int num_keypoints, + __global const float* Lx, + __global const float* Ly, + __global const float* Lt, + int rows, int cols, float octave_ratio, + int descriptor_size, + __global uchar* output_descriptors) +{ + int i = get_global_id(0); + + if (i >= num_keypoints) + return; + + float kpt_x = keypoints_x[i]; + float kpt_y = keypoints_y[i]; + float kpt_size = keypoints_size[i]; + + // cvRound implementation: rounds to nearest, half away from zero + float scale_f = 0.5f * kpt_size / octave_ratio; + int scale = (int)(scale_f + (scale_f >= 0 ? 0.5f : -0.5f)); + + float xf = kpt_x / octave_ratio; + float yf = kpt_y / octave_ratio; + + // Pattern size (default 10) + const int pattern_size = 10; + + // Sample steps for 3 grids: 2x2, 3x3, 4x4 + const int sample_step[3] = { + pattern_size, + (pattern_size * 2 + 2) / 3, // divUp(pattern_size * 2, 3) + pattern_size / 2 + }; + + // Buffer for M-LDB descriptor values (max 16 cells * 3 channels) + float values[48]; + + __global uchar* desc = &output_descriptors[i * descriptor_size]; + + // Initialize descriptor to zero + for (int b = 0; b < descriptor_size; b++) + desc[b] = 0; + + int dcount1 = 0; + + // For the three grids (2x2, 3x3, 4x4) + for (int z = 0; z < 3; z++) { + int dcount2 = 0; + const int step = sample_step[z]; + + for (int i = -pattern_size; i < pattern_size; i += step) { + for (int j = -pattern_size; j < pattern_size; j += step) { + float di = 0.0f, dx = 0.0f, dy = 0.0f; + int nsamples = 0; + + // Sample within each cell + for (int k = 0; k < step; k++) { + for (int l = 0; l < step; l++) { + // Get the coordinates of the sample point + const float sample_y = yf + (l + j) * scale; + const float sample_x = xf + (k + i) * scale; + + // cvRound implementation for sampling coordinates + const int y1 = (int)(sample_y + (sample_y >= 0 ? 0.5f : -0.5f)); + const int x1 = (int)(sample_x + (sample_x >= 0 ? 0.5f : -0.5f)); + + if (y1 < 0 || y1 >= rows || x1 < 0 || x1 >= cols) + continue; // Boundaries + + const int idx = y1 * cols + x1; + + const float ri = Lt[idx]; + const float rx = Lx[idx]; + const float ry = Ly[idx]; + + di += ri; + dx += rx; + dy += ry; + nsamples++; + } + } + + if (nsamples > 0) { + const float nsamples_inv = 1.0f / nsamples; + di *= nsamples_inv; + dx *= nsamples_inv; + dy *= nsamples_inv; + } + + // Store values (3 channels: Lt, Lx, Ly) + values[dcount2 * 3] = di; + values[dcount2 * 3 + 1] = dx; + values[dcount2 * 3 + 2] = dy; + dcount2++; + } + } + + // Do binary comparison for this grid + const int num = (z + 2) * (z + 2); + const int chan = 3; + + // Apply CV_TOGGLE_FLT to handle signed floats correctly + // This toggles the sign bit to allow correct integer comparison of floats + int* ivalues = (int*)values; + for (int i = 0; i < num * chan; i++) { + ivalues[i] = ivalues[i] ^ (ivalues[i] < 0 ? 0x7fffffff : 0); + } + + // Match CPU comparison order: iterate Cell FIRST, then Channel + // This produces: [Cell0-Ch0, Cell0-Ch1, Cell0-Ch2, Cell1-Ch0, Cell1-Ch1, ...] + for (int i = 0; i < num; i++) { + for (int j = i + 1; j < num; j++) { + for (int pos = 0; pos < chan; pos++) { + int ival = ivalues[chan * i + pos]; + if (ival > ivalues[chan * j + pos]) { + desc[dcount1 / 8] |= (1 << (dcount1 % 8)); + } + dcount1++; + } + } + } + } +} diff --git a/modules/features2d/src/opencl/kaze.cl b/modules/features2d/src/opencl/kaze.cl new file mode 100644 index 0000000000..4bc6169651 --- /dev/null +++ b/modules/features2d/src/opencl/kaze.cl @@ -0,0 +1,296 @@ +// 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 + +inline float gaussian(float x, float y, float sigma) +{ + return exp(-(x*x + y*y) / (2.0f*sigma*sigma)); +} + +/** + * @brief Compute KAZE upright 64-dimensional descriptor + * @details Matches CPU Get_KAZE_Upright_Descriptor_64 exactly: + * - integer scale via round(kpt_size/2) + * - bilinear with y1=(int)(y-0.5), y2=(int)(y+0.5) convention + * - clamp-on-border (not skip) + * + * @param Lx First-order x-derivative (sigma_size-scaled), row-major float array + * @param Ly First-order y-derivative (sigma_size-scaled), row-major float array + * @param lx_step Row stride of Lx/Ly in float elements (= cols for continuous mat) + * @param lx_rows Image height + * @param lx_cols Image width + * @param keypoints_x Keypoint x coordinates + * @param keypoints_y Keypoint y coordinates + * @param keypoints_size Keypoint sizes (diameter = 2*sigma) + * @param descriptors Output descriptor matrix, shape [nkeypoints, 64] + * @param nkeypoints Number of keypoints + */ +__kernel void +KAZE_compute_upright_descriptor_64( + __global const float* Lx, + __global const float* Ly, + int lx_step, int lx_rows, int lx_cols, + __global const float* keypoints_x, + __global const float* keypoints_y, + __global const float* keypoints_size, + __global float* descriptors, + int nkeypoints) +{ + int idx = get_global_id(0); + if (idx >= nkeypoints) + return; + + const int dsize = 64; + const int sample_step = 5; + const int pattern_size = 12; + + float kpt_x = keypoints_x[idx]; + float kpt_y = keypoints_y[idx]; + float kpt_sz = keypoints_size[idx]; + + int scale = (int)(kpt_sz / 2.0f + 0.5f); // cvRound equivalent + float xf = kpt_x; + float yf = kpt_y; + + float dx = 0.0f, dy = 0.0f, mdx = 0.0f, mdy = 0.0f; + float gauss_s1, gauss_s2; + float rx, ry; + float sample_x, sample_y; + int x1, y1, x2, y2; + int kx, ky, i, j, dcount = 0; + float fx, fy; + float res1, res2, res3, res4; + float len = 0.0f; + + float cx = -0.5f, cy = 0.5f; + + i = -8; + while (i < pattern_size) + { + j = -8; + i = i - 4; + cx += 1.0f; + cy = -0.5f; + + while (j < pattern_size) + { + dx = dy = mdx = mdy = 0.0f; + cy += 1.0f; + j = j - 4; + + ky = i + sample_step; + kx = j + sample_step; + + float ys = yf + (float)(ky * scale); + float xs = xf + (float)(kx * scale); + + for (int k = i; k < i + 9; k++) + { + for (int l = j; l < j + 9; l++) + { + sample_y = (float)k * scale + yf; + sample_x = (float)l * scale + xf; + + gauss_s1 = gaussian(xs - sample_x, ys - sample_y, 2.5f * scale); + + // Match CPU bilinear convention: y1=(int)(y-0.5), y2=(int)(y+0.5) + y1 = (int)(sample_y - 0.5f); + x1 = (int)(sample_x - 0.5f); + y1 = clamp(y1, 0, lx_rows - 1); + x1 = clamp(x1, 0, lx_cols - 1); + + y2 = (int)(sample_y + 0.5f); + x2 = (int)(sample_x + 0.5f); + y2 = clamp(y2, 0, lx_rows - 1); + x2 = clamp(x2, 0, lx_cols - 1); + + fx = sample_x - (float)x1; + fy = sample_y - (float)y1; + + res1 = Lx[y1 * lx_step + x1]; + res2 = Lx[y1 * lx_step + x2]; + res3 = Lx[y2 * lx_step + x1]; + res4 = Lx[y2 * lx_step + x2]; + rx = (1.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + + (1.0f - fx)*fy*res3 + fx*fy*res4; + + res1 = Ly[y1 * lx_step + x1]; + res2 = Ly[y1 * lx_step + x2]; + res3 = Ly[y2 * lx_step + x1]; + res4 = Ly[y2 * lx_step + x2]; + ry = (1.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + + (1.0f - fx)*fy*res3 + fx*fy*res4; + + rx = gauss_s1 * rx; + ry = gauss_s1 * ry; + + dx += rx; + dy += ry; + mdx += fabs(rx); + mdy += fabs(ry); + } + } + + gauss_s2 = gaussian(cx - 2.0f, cy - 2.0f, 1.5f); + + descriptors[idx * dsize + dcount++] = dx * gauss_s2; + descriptors[idx * dsize + dcount++] = dy * gauss_s2; + descriptors[idx * dsize + dcount++] = mdx * gauss_s2; + descriptors[idx * dsize + dcount++] = mdy * gauss_s2; + + len += (dx*dx + dy*dy + mdx*mdx + mdy*mdy) * gauss_s2 * gauss_s2; + + j += 9; + } + i += 9; + } + + // L2 normalize + len = sqrt(len); + if (len > 1e-10f) + { + float len_inv = 1.0f / len; + for (i = 0; i < dsize; i++) + descriptors[idx * dsize + i] *= len_inv; + } +} + +/** + * @brief Compute KAZE upright 128-dimensional descriptor (extended mode) + * @details Matches CPU Get_KAZE_Upright_Descriptor_128: splits dx/dy by sign of the + * cross-component for 8 values per subregion. + */ +__kernel void +KAZE_compute_upright_descriptor_128( + __global const float* Lx, + __global const float* Ly, + int lx_step, int lx_rows, int lx_cols, + __global const float* keypoints_x, + __global const float* keypoints_y, + __global const float* keypoints_size, + __global float* descriptors, + int nkeypoints) +{ + int idx = get_global_id(0); + if (idx >= nkeypoints) + return; + + const int dsize = 128; + const int sample_step = 5; + const int pattern_size = 12; + + float kpt_x = keypoints_x[idx]; + float kpt_y = keypoints_y[idx]; + float kpt_sz = keypoints_size[idx]; + + int scale = (int)(kpt_sz / 2.0f + 0.5f); + float xf = kpt_x; + float yf = kpt_y; + + float dxp = 0.0f, dxn = 0.0f, mdxp = 0.0f, mdxn = 0.0f; + float dyp = 0.0f, dyn = 0.0f, mdyp = 0.0f, mdyn = 0.0f; + float gauss_s1, gauss_s2; + float rx, ry; + float sample_x, sample_y; + int x1, y1, x2, y2; + int kx, ky, i, j, dcount = 0; + float fx, fy; + float res1, res2, res3, res4; + float len = 0.0f; + + float cx = -0.5f, cy = 0.5f; + + i = -8; + while (i < pattern_size) + { + j = -8; + i = i - 4; + cx += 1.0f; + cy = -0.5f; + + while (j < pattern_size) + { + dxp = dxn = mdxp = mdxn = 0.0f; + dyp = dyn = mdyp = mdyn = 0.0f; + cy += 1.0f; + j = j - 4; + + ky = i + sample_step; + kx = j + sample_step; + + float ys = yf + (float)(ky * scale); + float xs = xf + (float)(kx * scale); + + for (int k = i; k < i + 9; k++) + { + for (int l = j; l < j + 9; l++) + { + sample_y = (float)k * scale + yf; + sample_x = (float)l * scale + xf; + + gauss_s1 = gaussian(xs - sample_x, ys - sample_y, 2.5f * scale); + + y1 = (int)(sample_y - 0.5f); + x1 = (int)(sample_x - 0.5f); + y1 = clamp(y1, 0, lx_rows - 1); + x1 = clamp(x1, 0, lx_cols - 1); + + y2 = (int)(sample_y + 0.5f); + x2 = (int)(sample_x + 0.5f); + y2 = clamp(y2, 0, lx_rows - 1); + x2 = clamp(x2, 0, lx_cols - 1); + + fx = sample_x - (float)x1; + fy = sample_y - (float)y1; + + res1 = Lx[y1 * lx_step + x1]; + res2 = Lx[y1 * lx_step + x2]; + res3 = Lx[y2 * lx_step + x1]; + res4 = Lx[y2 * lx_step + x2]; + rx = (1.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + + (1.0f - fx)*fy*res3 + fx*fy*res4; + + res1 = Ly[y1 * lx_step + x1]; + res2 = Ly[y1 * lx_step + x2]; + res3 = Ly[y2 * lx_step + x1]; + res4 = Ly[y2 * lx_step + x2]; + ry = (1.0f - fx)*(1.0f - fy)*res1 + fx*(1.0f - fy)*res2 + + (1.0f - fx)*fy*res3 + fx*fy*res4; + + rx = gauss_s1 * rx; + ry = gauss_s1 * ry; + + if (ry >= 0.0f) { dxp += rx; mdxp += fabs(rx); } + else { dxn += rx; mdxn += fabs(rx); } + if (rx >= 0.0f) { dyp += ry; mdyp += fabs(ry); } + else { dyn += ry; mdyn += fabs(ry); } + } + } + + gauss_s2 = gaussian(cx - 2.0f, cy - 2.0f, 1.5f); + + descriptors[idx * dsize + dcount++] = dxp * gauss_s2; + descriptors[idx * dsize + dcount++] = dxn * gauss_s2; + descriptors[idx * dsize + dcount++] = mdxp * gauss_s2; + descriptors[idx * dsize + dcount++] = mdxn * gauss_s2; + descriptors[idx * dsize + dcount++] = dyp * gauss_s2; + descriptors[idx * dsize + dcount++] = dyn * gauss_s2; + descriptors[idx * dsize + dcount++] = mdyp * gauss_s2; + descriptors[idx * dsize + dcount++] = mdyn * gauss_s2; + + len += (dxp*dxp + dxn*dxn + mdxp*mdxp + mdxn*mdxn + + dyp*dyp + dyn*dyn + mdyp*mdyp + mdyn*mdyn) * gauss_s2 * gauss_s2; + + j += 9; + } + i += 9; + } + + len = sqrt(len); + if (len > 1e-10f) + { + float len_inv = 1.0f / len; + for (i = 0; i < dsize; i++) + descriptors[idx * dsize + i] *= len_inv; + } +} diff --git a/modules/features2d/test/ocl/test_akaze.cpp b/modules/features2d/test/ocl/test_akaze.cpp new file mode 100644 index 0000000000..55ae155049 --- /dev/null +++ b/modules/features2d/test/ocl/test_akaze.cpp @@ -0,0 +1,140 @@ +// 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" +#include "cvconfig.h" +#include "opencv2/ts/ocl_test.hpp" + +#ifdef HAVE_OPENCL + +namespace opencv_test { +namespace ocl { +TEST(Features2d_AKAZE, ocl_accuracy) +{ + Mat testImg(640, 480, CV_8U); + theRNG().fill(testImg, RNG::UNIFORM, Scalar(0), Scalar(255), true); + + // CPU version - use MLDB_UPRIGHT to match GPU implementation + Ptr akaze_cpu = AKAZE::create(AKAZE::DESCRIPTOR_MLDB_UPRIGHT, 0, 3, 0.001f, 1, 1, KAZE::DIFF_PM_G2); + vector kp_cpu; + Mat desc_cpu; + akaze_cpu->detectAndCompute(testImg, noArray(), kp_cpu, desc_cpu); + + // OpenCL version - use MLDB_UPRIGHT to match GPU implementation + UMat testImg_umat; + testImg.copyTo(testImg_umat); + Ptr akaze_ocl = AKAZE::create(AKAZE::DESCRIPTOR_MLDB_UPRIGHT, 0, 3, 0.001f, 1, 1, KAZE::DIFF_PM_G2); + vector kp_ocl; + UMat desc_ocl_umat; + akaze_ocl->detectAndCompute(testImg_umat, noArray(), kp_ocl, desc_ocl_umat); + Mat desc_ocl = desc_ocl_umat.getMat(ACCESS_READ); + + // Check that both detected keypoints + ASSERT_FALSE(kp_cpu.empty()) << "CPU should detect keypoints"; + ASSERT_FALSE(kp_ocl.empty()) << "OpenCL should detect keypoints"; + + // Allow small keypoint count difference due to border handling + float kp_ratio = (float)kp_ocl.size() / kp_cpu.size(); + EXPECT_GT(kp_ratio, 0.95f) << "OpenCL keypoint count should be within 5% of CPU, got " << kp_ratio; + EXPECT_LT(kp_ratio, 1.05f) << "OpenCL keypoint count should be within 5% of CPU, got " << kp_ratio; + + // Check descriptor dimensions match + ASSERT_EQ(desc_cpu.cols, desc_ocl.cols) << "Descriptor size should match"; + ASSERT_EQ(desc_cpu.type(), desc_ocl.type()) << "Descriptor type should match"; + + // Match keypoints by position (within 1 pixel tolerance) + int matched_kpts = 0; + int total_compared = 0; + int matching_bytes = 0; + + for (size_t i = 0; i < kp_cpu.size(); i++) { + // Find matching keypoint in OpenCL results + for (size_t j = 0; j < kp_ocl.size(); j++) { + float dx = fabs(kp_cpu[i].pt.x - kp_ocl[j].pt.x); + float dy = fabs(kp_cpu[i].pt.y - kp_ocl[j].pt.y); + + if (dx < 1.0f && dy < 1.0f) { + // Found matching keypoint, compare descriptors + matched_kpts++; + for (int k = 0; k < desc_cpu.cols; k++) { + if (desc_cpu.at((int)i, k) == desc_ocl.at((int)j, k)) { + matching_bytes++; + } + } + total_compared += desc_cpu.cols; + break; + } + } + } + + float match_rate = total_compared > 0 ? (float)matching_bytes / total_compared : 0.0f; + EXPECT_GT(matched_kpts, (int)(kp_cpu.size() * 0.9f)) << "Should match at least 90% of keypoints by position"; + EXPECT_GT(match_rate, 0.95f) << "Descriptor match rate should be > 95%, got " << match_rate; +} + +TEST(Features2d_KAZE, ocl_accuracy) +{ + Mat testImg(640, 480, CV_8U); + theRNG().fill(testImg, RNG::UNIFORM, Scalar(0), Scalar(255), true); + + // CPU version - use upright=true to match GPU implementation + Ptr kaze_cpu = KAZE::create(false, true, 0.001f, 4, 4, KAZE::DIFF_PM_G2); + vector kp_cpu; + Mat desc_cpu; + kaze_cpu->detectAndCompute(testImg, noArray(), kp_cpu, desc_cpu); + + // OpenCL version - use upright=true to match GPU implementation + UMat testImg_umat; + testImg.copyTo(testImg_umat); + Ptr kaze_ocl = KAZE::create(false, true, 0.001f, 4, 4, KAZE::DIFF_PM_G2); + vector kp_ocl; + UMat desc_ocl_umat; + kaze_ocl->detectAndCompute(testImg_umat, noArray(), kp_ocl, desc_ocl_umat); + Mat desc_ocl = desc_ocl_umat.getMat(ACCESS_READ); + + // Check that both detected keypoints + ASSERT_FALSE(kp_cpu.empty()) << "CPU should detect keypoints"; + ASSERT_FALSE(kp_ocl.empty()) << "OpenCL should detect keypoints"; + + // Allow small keypoint count difference due to border handling + float kp_ratio = (float)kp_ocl.size() / kp_cpu.size(); + EXPECT_GT(kp_ratio, 0.95f) << "OpenCL keypoint count should be within 5% of CPU, got " << kp_ratio; + EXPECT_LT(kp_ratio, 1.05f) << "OpenCL keypoint count should be within 5% of CPU, got " << kp_ratio; + + // Check descriptor dimensions match + ASSERT_EQ(desc_cpu.cols, desc_ocl.cols) << "Descriptor size should match"; + ASSERT_EQ(desc_cpu.type(), desc_ocl.type()) << "Descriptor type should match"; + + // Match keypoints by position (within 1 pixel tolerance) + int matched_kpts = 0; + int total_compared = 0; + double total_diff = 0; + + for (size_t i = 0; i < kp_cpu.size(); i++) { + // Find matching keypoint in OpenCL results + for (size_t j = 0; j < kp_ocl.size(); j++) { + float dx = fabs(kp_cpu[i].pt.x - kp_ocl[j].pt.x); + float dy = fabs(kp_cpu[i].pt.y - kp_ocl[j].pt.y); + + if (dx < 1.0f && dy < 1.0f) { + // Found matching keypoint, compare descriptors + matched_kpts++; + for (int k = 0; k < desc_cpu.cols; k++) { + double diff = fabs(desc_cpu.at((int)i, k) - desc_ocl.at((int)j, k)); + total_diff += diff; + total_compared++; + } + break; + } + } + } + + float avg_diff = total_compared > 0 ? (float)(total_diff / total_compared) : 0.0f; + EXPECT_GT(matched_kpts, (int)(kp_cpu.size() * 0.9f)) << "Should match at least 90% of keypoints by position"; + EXPECT_LT(avg_diff, 0.01f) << "Average descriptor difference should be < 0.01, got " << avg_diff; +} + +}} // namespace + +#endif diff --git a/modules/features2d/test/test_akaze.cpp b/modules/features2d/test/test_akaze.cpp index e24894129a..c4234ab28f 100644 --- a/modules/features2d/test/test_akaze.cpp +++ b/modules/features2d/test/test_akaze.cpp @@ -9,8 +9,7 @@ namespace opencv_test { namespace { 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); + theRNG().fill(testImg, RNG::UNIFORM, Scalar(0), Scalar(255), true); Ptr ext = AKAZE::create(AKAZE::DESCRIPTOR_MLDB, 0, 3, 0.001f, 1, 1, KAZE::DIFF_PM_G2); vector detAndCompKps; @@ -49,8 +48,7 @@ TEST(Features2d_AKAZE, uninitialized_and_nans) TEST(Features2d_KAZE, diffusivity_charbonnier) { Mat testImg(200, 200, CV_8U); - RNG rng(42); - rng.fill(testImg, RNG::UNIFORM, Scalar(0), Scalar(255), true); + theRNG().fill(testImg, RNG::UNIFORM, Scalar(0), Scalar(255), true); // KAZE with DIFF_CHARBONNIER Ptr kaze_charbonnier = KAZE::create(false, false, 0.001f, 4, 4, KAZE::DIFF_CHARBONNIER);