From 4ebaf2dbb4ad1c1420a1c69858441540654e7ab0 Mon Sep 17 00:00:00 2001 From: gideok Kim Date: Tue, 10 Mar 2026 23:18:02 +0900 Subject: [PATCH 01/19] imgproc: fix fitEllipseDirect/AMS determinant threshold for near-circular data MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit fitEllipseDirect: replace det(M) threshold with eigenvector quality check - Add Ts ≈ 0 guard to avoid division by zero in Schur complement - Move eigenNonSymmetric inside perturbation loop - Validate eigenvector with 4ac-b² > 1e-6*||v||² to filter garbage from complex eigenvalues fitEllipseAMS: scale threshold by 1/n^5 to match det(M) magnitude Add Imgproc_FitEllipseDirect_NearCircular regression test --- modules/imgproc/src/shapedescr.cpp | 32 +++++++++++++------ modules/imgproc/test/test_fitellipse_ams.cpp | 20 ++++++++++++ .../imgproc/test/test_fitellipse_direct.cpp | 22 +++++++++++++ 3 files changed, 64 insertions(+), 10 deletions(-) diff --git a/modules/imgproc/src/shapedescr.cpp b/modules/imgproc/src/shapedescr.cpp index 6ddd855e32..4952a9722e 100644 --- a/modules/imgproc/src/shapedescr.cpp +++ b/modules/imgproc/src/shapedescr.cpp @@ -579,7 +579,9 @@ cv::RotatedRect cv::fitEllipseAMS( InputArray _points ) M(4,3)=DM(3,4); M(4,4)=DM(4,4); - if (fabs(cv::determinant(M)) > 1.0e-10) { + double npow = (double)n * (double)n; + npow = npow * npow * (double)n; // n^5 + if (fabs(cv::determinant(M)) > 1.0e-10 / npow) { break; } @@ -703,6 +705,8 @@ cv::RotatedRect cv::fitEllipseDirect( InputArray _points ) Matx pVec; double x0, y0, a, b, theta, Ts; + Mat eVal, eVec; + double cond[3]; double s = 0; for( i = 0; i < n; i++ ) @@ -763,6 +767,11 @@ cv::RotatedRect cv::fitEllipseDirect( InputArray _points ) Ts=(-(DM(3,5)*DM(4,4)*DM(5,3)) + DM(3,4)*DM(4,5)*DM(5,3) + DM(3,5)*DM(4,3)*DM(5,4) - \ DM(3,3)*DM(4,5)*DM(5,4) - DM(3,4)*DM(4,3)*DM(5,5) + DM(3,3)*DM(4,4)*DM(5,5)); + if (fabs(Ts) < DBL_EPSILON) { + eps = (float)(s/(n*2)*1e-2); + continue; + } + M(0,0) = (DM(2,0) + (DM(2,3)*TM(0,0) + DM(2,4)*TM(1,0) + DM(2,5)*TM(2,0))/Ts)/2.; M(0,1) = (DM(2,1) + (DM(2,3)*TM(0,1) + DM(2,4)*TM(1,1) + DM(2,5)*TM(2,1))/Ts)/2.; M(0,2) = (DM(2,2) + (DM(2,3)*TM(0,2) + DM(2,4)*TM(1,2) + DM(2,5)*TM(2,2))/Ts)/2.; @@ -773,18 +782,9 @@ cv::RotatedRect cv::fitEllipseDirect( InputArray _points ) M(2,1) = (DM(0,1) + (DM(0,3)*TM(0,1) + DM(0,4)*TM(1,1) + DM(0,5)*TM(2,1))/Ts)/2.; M(2,2) = (DM(0,2) + (DM(0,3)*TM(0,2) + DM(0,4)*TM(1,2) + DM(0,5)*TM(2,2))/Ts)/2.; - double det = cv::determinant(M); - if (fabs(det) > 1.0e-10) - break; - eps = (float)(s/(n*2)*1e-2); - } - - if( iter < 2 ) { - Mat eVal, eVec; eigenNonSymmetric(M, eVal, eVec); // Select the eigen vector {a,b,c} which satisfies 4ac-b^2 > 0 - double cond[3]; cond[0]=(4.0 * eVec.at(0,0) * eVec.at(0,2) - eVec.at(0,1) * eVec.at(0,1)); cond[1]=(4.0 * eVec.at(1,0) * eVec.at(1,2) - eVec.at(1,1) * eVec.at(1,1)); cond[2]=(4.0 * eVec.at(2,0) * eVec.at(2,2) - eVec.at(2,1) * eVec.at(2,1)); @@ -793,6 +793,18 @@ cv::RotatedRect cv::fitEllipseDirect( InputArray _points ) } else { i = (cond[0](i,0), v1 = eVec.at(i,1), v2 = eVec.at(i,2); + double vnorm2 = v0*v0 + v1*v1 + v2*v2; + if (cond[i] > 1e-6 * vnorm2) + break; + } + eps = (float)(s/(n*2)*1e-2); + } + + if( iter < 2 ) { double norm = std::sqrt(eVec.at(i,0)*eVec.at(i,0) + eVec.at(i,1)*eVec.at(i,1) + eVec.at(i,2)*eVec.at(i,2)); if (((eVec.at(i,0)<0.0 ? -1 : 1) * (eVec.at(i,1)<0.0 ? -1 : 1) * (eVec.at(i,2)<0.0 ? -1 : 1)) <= 0.0) { norm=-1.0*norm; diff --git a/modules/imgproc/test/test_fitellipse_ams.cpp b/modules/imgproc/test/test_fitellipse_ams.cpp index f2c9d1f793..b88660e1f3 100644 --- a/modules/imgproc/test/test_fitellipse_ams.cpp +++ b/modules/imgproc/test/test_fitellipse_ams.cpp @@ -337,6 +337,26 @@ TEST(Imgproc_FitEllipseAMS_Issue_7, accuracy) { EXPECT_TRUE(checkEllipse(ellipseAMSTest, ellipseAMSTrue, tol)); } +TEST(Imgproc_FitEllipseAMS_NearCircular, accuracy) +{ + std::vector points; + double cx = 27.0, cy = 27.0, a = 17.0, b = 16.5; + for (int i = 0; i < 360; i++) { + double theta = 2.0 * CV_PI * i / 360.0; + points.push_back(cv::Point2f( + (float)(cx + a * cos(theta)), + (float)(cy + b * sin(theta)))); + } + + cv::RotatedRect ams = cv::fitEllipseAMS(points); + + // AMS should produce a valid result close to ground truth + EXPECT_NEAR(ams.center.x, 27.0, 0.5); + EXPECT_NEAR(ams.center.y, 27.0, 0.5); + EXPECT_NEAR(std::max(ams.size.width, ams.size.height), 34.0, 1.0); + EXPECT_NEAR(std::min(ams.size.width, ams.size.height), 33.0, 1.0); +} + TEST(Imgproc_FitEllipseAMS_HorizontalLine, accuracy) { vector pts({{-300, 100}, {-200, 100}, {-100, 100}, {0, 100}, {100, 100}, {200, 100}, {300, 100}}); const RotatedRect el = fitEllipseAMS(pts); diff --git a/modules/imgproc/test/test_fitellipse_direct.cpp b/modules/imgproc/test/test_fitellipse_direct.cpp index e41c52764d..81d05a9280 100644 --- a/modules/imgproc/test/test_fitellipse_direct.cpp +++ b/modules/imgproc/test/test_fitellipse_direct.cpp @@ -337,6 +337,28 @@ TEST(Imgproc_FitEllipseDirect_Issue_7, accuracy) { EXPECT_TRUE(checkEllipse(ellipseDirectTest, ellipseDirectTrue, tol)); } +TEST(Imgproc_FitEllipseDirect_NearCircular, accuracy) +{ + // 360 points on a near-circular ellipse (a=17, b=16.5) + // This data previously triggered unnecessary fallback to fitEllipseNoDirect + std::vector points; + double cx = 27.0, cy = 27.0, a = 17.0, b = 16.5; + for (int i = 0; i < 360; i++) { + double theta = 2.0 * CV_PI * i / 360.0; + points.push_back(cv::Point2f( + (float)(cx + a * cos(theta)), + (float)(cy + b * sin(theta)))); + } + + cv::RotatedRect direct = cv::fitEllipseDirect(points); + + // Direct should produce a valid result close to ground truth + EXPECT_NEAR(direct.center.x, 27.0, 0.1); + EXPECT_NEAR(direct.center.y, 27.0, 0.1); + EXPECT_NEAR(std::max(direct.size.width, direct.size.height), 34.0, 0.5); + EXPECT_NEAR(std::min(direct.size.width, direct.size.height), 33.0, 0.5); +} + TEST(Imgproc_FitEllipseDirect_HorizontalLine, accuracy) { vector pts({{-300, 100}, {-200, 100}, {-100, 100}, {0, 100}, {100, 100}, {200, 100}, {300, 100}}); const RotatedRect el = fitEllipseDirect(pts); From 4acd2ed6cc21bcd974a4e5183f7ce2f2abdc78e1 Mon Sep 17 00:00:00 2001 From: Anand Mahesh Date: Sun, 29 Mar 2026 21:43:23 +0530 Subject: [PATCH 02/19] features2d: add OpenCL acceleration for BFMatcher cross-check MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit BFMatcher::match() with crossCheck=true previously skipped the OCL dispatch in knnMatchImpl entirely, falling back to CPU even when UMat inputs and an OpenCL device were available. This adds ocl_matchWithCrossCheck() for CV_32FC1 descriptors (e.g. SIFT, SURF): both the forward and reverse nearest-neighbour passes run on the GPU via the existing ocl_matchSingle() kernel, then the cross-check filter runs on the CPU. Only two small index arrays (1×N ints) are downloaded — the O(N²×D) distance work stays on the device. The OCL dispatch in knnMatchImpl is also refactored to unify the Mat/UMat train collection selection before branching on crossCheck. On a NVIDIA RTX 3060 with SIFT descriptors the OCL path is 6–9× faster than CPU at 2k–10k features per image. Binary descriptors (ORB, BRIEF — CV_8U) are unaffected; the existing type guard in ocl_matchSingle keeps them on the CPU path as before. Also adds a correctness test (Features2d_BFMatcher_CrossCheck) and an OCL perf test (BruteForceMatcherFixture/MatchCrossCheck). --- .../perf/opencl/perf_brute_force_matcher.cpp | 21 ++++++ modules/features2d/src/matchers.cpp | 65 +++++++++++++++---- .../test/test_matchers_algorithmic.cpp | 43 ++++++++++++ 3 files changed, 115 insertions(+), 14 deletions(-) diff --git a/modules/features2d/perf/opencl/perf_brute_force_matcher.cpp b/modules/features2d/perf/opencl/perf_brute_force_matcher.cpp index 439140a5ff..6af06d37ad 100644 --- a/modules/features2d/perf/opencl/perf_brute_force_matcher.cpp +++ b/modules/features2d/perf/opencl/perf_brute_force_matcher.cpp @@ -123,6 +123,27 @@ OCL_PERF_TEST_P(BruteForceMatcherFixture, RadiusMatch, ::testing::Combine(OCL_PE SANITY_CHECK_MATCHES(matches1, 1e-3); } +OCL_PERF_TEST_P(BruteForceMatcherFixture, MatchCrossCheck, ::testing::Combine(OCL_PERF_ENUM(OCL_SIZE_1, OCL_SIZE_2, OCL_SIZE_3), OCL_PERF_ENUM((MatType)CV_32FC1) ) ) +{ + const Size_MatType_t params = GetParam(); + const Size srcSize = get<0>(params); + const int type = get<1>(params); + + checkDeviceMaxMemoryAllocSize(srcSize, type); + + vector matches; + UMat uquery(srcSize, type), utrain(srcSize, type); + + declare.in(uquery, utrain, WARMUP_RNG); + + BFMatcher matcher(NORM_L2, true /*crossCheck*/); + + OCL_TEST_CYCLE() + matcher.match(uquery, utrain, matches); + + SANITY_CHECK_MATCHES(matches, 1e-3); +} + } // ocl } // cvtest diff --git a/modules/features2d/src/matchers.cpp b/modules/features2d/src/matchers.cpp index 97d0a02417..950fd96d58 100644 --- a/modules/features2d/src/matchers.cpp +++ b/modules/features2d/src/matchers.cpp @@ -752,6 +752,49 @@ static bool ocl_knnMatch(InputArray query, InputArray _train, std::vector< std:: return false; return true; } + +// Run match in both directions on GPU, cross-check filter on CPU. +// Reuses ocl_matchSingle so no new kernel is needed. +static bool ocl_matchWithCrossCheck(InputArray query, InputArray train, + std::vector< std::vector >& matches, int dstType) +{ + // Forward pass: for each query descriptor find nearest in train + UMat fwdIdx, fwdDist; + if (!ocl_matchSingle(query, train, fwdIdx, fwdDist, dstType)) + return false; + + // Reverse pass: for each train descriptor find nearest in query + UMat revIdx, revDist; + if (!ocl_matchSingle(train, query, revIdx, revDist, dstType)) + return false; + + // Download index arrays (1 x N ints each — cheap) + Mat fwdIdxCPU = fwdIdx.getMat(ACCESS_READ); + Mat revIdxCPU = revIdx.getMat(ACCESS_READ); + Mat fwdDistCPU = fwdDist.getMat(ACCESS_READ); + + if (fwdIdxCPU.empty() || revIdxCPU.empty() || fwdDistCPU.empty()) + return false; + + const int nQuery = fwdIdxCPU.cols; + const int nTrain = revIdxCPU.cols; + const int* fwd = fwdIdxCPU.ptr(); + const int* rev = revIdxCPU.ptr(); + const float* dist = fwdDistCPU.ptr(); + + matches.clear(); + matches.reserve(nQuery); + + for (int q = 0; q < nQuery; ++q) + { + int t = fwd[q]; + if (t >= 0 && t < nTrain && rev[t] == q) + { + matches.push_back(std::vector(1, DMatch(q, t, 0, dist[q]))); + } + } + return true; +} #endif void BFMatcher::knnMatchImpl( InputArray _queryDescriptors, std::vector >& matches, int knn, @@ -794,21 +837,15 @@ void BFMatcher::knnMatchImpl( InputArray _queryDescriptors, std::vectorknnMatch(usources, utargets, match, 1, mask, true)); } +// Verify that cross-check BFMatcher gives identical results via Mat (CPU) and UMat (OCL or CPU). +// When OpenCL is active the UMat path exercises ocl_matchWithCrossCheck; when it is not, +// both paths fall through to the same CPU code — either way the results must match. +TEST(Features2d_BFMatcher_CrossCheck, ocl_matches_cpu) +{ + RNG rng(42); + const int nQuery = 200; + const int nTrain = 400; + const int dim = 128; + + // Float descriptors: the OCL dispatch in knnMatchImpl requires CV_32FC1 + Mat queryMat(nQuery, dim, CV_32FC1); + Mat trainMat(nTrain, dim, CV_32FC1); + rng.fill(queryMat, RNG::UNIFORM, 0.f, 1.f); + rng.fill(trainMat, RNG::UNIFORM, 0.f, 1.f); + + // CPU reference: Mat inputs always take the CPU path + Ptr cpuMatcher = BFMatcher::create(NORM_L2, true /*crossCheck*/); + vector cpuMatches; + cpuMatcher->match(queryMat, trainMat, cpuMatches); + + // UMat path: activates OCL dispatch when OpenCL is available + UMat queryUMat = queryMat.getUMat(ACCESS_READ); + UMat trainUMat = trainMat.getUMat(ACCESS_READ); + Ptr oclMatcher = BFMatcher::create(NORM_L2, true /*crossCheck*/); + vector oclMatches; + oclMatcher->match(queryUMat, trainUMat, oclMatches); + + // Both paths must return the same set of matches (order may differ) + ASSERT_EQ(cpuMatches.size(), oclMatches.size()); + + auto byQuery = [](const DMatch& a, const DMatch& b) { return a.queryIdx < b.queryIdx; }; + sort(cpuMatches.begin(), cpuMatches.end(), byQuery); + sort(oclMatches.begin(), oclMatches.end(), byQuery); + + for (size_t i = 0; i < cpuMatches.size(); ++i) + { + EXPECT_EQ(cpuMatches[i].queryIdx, oclMatches[i].queryIdx) << "at index " << i; + EXPECT_EQ(cpuMatches[i].trainIdx, oclMatches[i].trainIdx) << "at index " << i; + EXPECT_NEAR(cpuMatches[i].distance, oclMatches[i].distance, 1e-3f) << "at index " << i; + } +} + }} // namespace From 758e8620accbe7fcd4f460e8fd73c584bd72542b Mon Sep 17 00:00:00 2001 From: Ismail Date: Sat, 4 Apr 2026 00:02:52 +0200 Subject: [PATCH 03/19] photo: remove redundant code in illuminationChange Fixes redundant code in Cloning::illuminationChange() which performed unnecessary copyTo operations. This satisfies Issue #22056. Testing: No functional changes made; logic identical to before. Resolves #22056 --- modules/photo/src/seamless_cloning_impl.cpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/modules/photo/src/seamless_cloning_impl.cpp b/modules/photo/src/seamless_cloning_impl.cpp index 377bdb8f8d..c2e5c13478 100644 --- a/modules/photo/src/seamless_cloning_impl.cpp +++ b/modules/photo/src/seamless_cloning_impl.cpp @@ -428,11 +428,6 @@ void Cloning::illuminationChange(Mat &I, Mat &mask, Mat &wmask, Mat &cloned, flo multiply(multY,multy_temp,patchGradientY); patchNaNs(patchGradientY); - Mat zeroMask = (patchGradientX != 0); - - patchGradientX.copyTo(patchGradientX, zeroMask); - patchGradientY.copyTo(patchGradientY, zeroMask); - evaluate(I,wmask,cloned); } From 1bef8b81f441e9a8f860fc965e5e937b68054954 Mon Sep 17 00:00:00 2001 From: Samaresh Kumar Singh Date: Tue, 7 Apr 2026 21:24:04 -0500 Subject: [PATCH 04/19] cmake: find and link libcblas separately for LAPACK/Generic On systems like OpenBSD, CBLAS functions are provided by a standalone libcblas that is not returned by find_package(LAPACK). This caused link failures when building libopencv_core because cblas_sgemm and friends from hal_internal.cpp could not be resolved. Fixes #28768 --- cmake/OpenCVFindLAPACK.cmake | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/cmake/OpenCVFindLAPACK.cmake b/cmake/OpenCVFindLAPACK.cmake index 3b1c435bdf..61a149c234 100644 --- a/cmake/OpenCVFindLAPACK.cmake +++ b/cmake/OpenCVFindLAPACK.cmake @@ -230,11 +230,18 @@ if(WITH_LAPACK) find_path(CBLAS_INCLUDE_DIR "cblas.h") endif() if(CBLAS_INCLUDE_DIR AND LAPACKE_INCLUDE_DIR) + if(NOT DEFINED CBLAS_LIBRARY) + find_library(CBLAS_LIBRARY cblas) + endif() + set(_lapack_generic_libs "${LAPACK_LIBRARIES}") + if(CBLAS_LIBRARY) + list(APPEND _lapack_generic_libs "${CBLAS_LIBRARY}") + endif() ocv_lapack_check(IMPL "LAPACK/Generic" CBLAS_H "cblas.h" LAPACKE_H "lapacke.h" INCLUDE_DIR "${CBLAS_INCLUDE_DIR}" "${LAPACKE_INCLUDE_DIR}" - LIBRARIES "${LAPACK_LIBRARIES}") + LIBRARIES "${_lapack_generic_libs}") elseif(APPLE) ocv_lapack_check(IMPL "LAPACK/Apple" CBLAS_H "Accelerate/Accelerate.h" From 3282bfc1498e1b8bc99664ef71b568f737c77ca2 Mon Sep 17 00:00:00 2001 From: Alexander Smorkalov Date: Thu, 9 Apr 2026 14:54:48 +0300 Subject: [PATCH 05/19] Added missing includes for CV_LOG_xxx. --- modules/calib3d/src/precomp.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/modules/calib3d/src/precomp.hpp b/modules/calib3d/src/precomp.hpp index 118e9dc974..e6b11705bb 100644 --- a/modules/calib3d/src/precomp.hpp +++ b/modules/calib3d/src/precomp.hpp @@ -42,6 +42,7 @@ #ifndef __OPENCV_PRECOMP_H__ #define __OPENCV_PRECOMP_H__ +#include #include "opencv2/core/utility.hpp" #include "opencv2/core/private.hpp" From 570e32cc9b004044ea675ca1ea8f156b1564b86d Mon Sep 17 00:00:00 2001 From: Andrew Yooeun Chun Date: Thu, 9 Apr 2026 22:03:14 +0900 Subject: [PATCH 06/19] videoio: fix FFmpeg VideoCapture rejecting CAP_PROP_FORMAT=CV_8UC3 --- modules/videoio/src/cap_ffmpeg_impl.hpp | 2 +- modules/videoio/test/test_ffmpeg.cpp | 16 ++++++++++++++++ 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/modules/videoio/src/cap_ffmpeg_impl.hpp b/modules/videoio/src/cap_ffmpeg_impl.hpp index 33f6eecc70..b892bd561b 100644 --- a/modules/videoio/src/cap_ffmpeg_impl.hpp +++ b/modules/videoio/src/cap_ffmpeg_impl.hpp @@ -1105,7 +1105,7 @@ bool CvCapture_FFMPEG::open(const char* _filename, int index, const Ptr Date: Fri, 10 Apr 2026 14:16:29 +0200 Subject: [PATCH 07/19] Merge pull request #28785 from chacha21:tiff_32F_compression Allow TIFF compression schemes for 32F#28785 Related to [https://github.com/opencv/opencv/issues/28775](https://github.com/opencv/opencv/issues/28775) Previously, only 32FC3+SGILOG could be specified for TIFF encoding. But when floats use quantization, compression schemes can be efficient even on 32F data. This PR will allow them. ### 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 - [X] There is a reference to the original bug report and related work - [ ] 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 I just don't know what kind of accuracy/performance tests should be added. --- .../imgcodecs/include/opencv2/imgcodecs.hpp | 9 +- modules/imgcodecs/src/grfmt_tiff.cpp | 83 ++++++++++++++++--- modules/imgcodecs/src/grfmt_tiff.hpp | 1 + modules/imgcodecs/test/test_tiff.cpp | 43 +++++++++- 4 files changed, 121 insertions(+), 15 deletions(-) diff --git a/modules/imgcodecs/include/opencv2/imgcodecs.hpp b/modules/imgcodecs/include/opencv2/imgcodecs.hpp index 6c59d55ec3..bb6fbeb998 100644 --- a/modules/imgcodecs/include/opencv2/imgcodecs.hpp +++ b/modules/imgcodecs/include/opencv2/imgcodecs.hpp @@ -109,7 +109,7 @@ enum ImwriteFlags { IMWRITE_TIFF_RESUNIT = 256,//!< For TIFF, use to specify which DPI resolution unit to set. See ImwriteTiffResolutionUnitFlags. Default is IMWRITE_TIFF_RESOLUTION_UNIT_INCH. IMWRITE_TIFF_XDPI = 257,//!< For TIFF, use to specify the X direction DPI IMWRITE_TIFF_YDPI = 258,//!< For TIFF, use to specify the Y direction DPI - IMWRITE_TIFF_COMPRESSION = 259,//!< For TIFF, use to specify the image compression scheme. See cv::ImwriteTiffCompressionFlags. Note, for images whose depth is CV_32F, only libtiff's SGILOG compression scheme is used. For other supported depths, the compression scheme can be specified by this flag; LZW compression is the default. + IMWRITE_TIFF_COMPRESSION = 259,//!< For TIFF, use to specify the image compression scheme. See cv::ImwriteTiffCompressionFlags. The compression scheme can be specified by this flag; the default is LZW compression, except for 32F depth where it is NONE IMWRITE_TIFF_ROWSPERSTRIP = 278,//!< For TIFF, use to specify the number of rows per strip. IMWRITE_TIFF_PREDICTOR = 317,//!< For TIFF, use to specify predictor. See cv::ImwriteTiffPredictorFlags. Default is IMWRITE_TIFF_PREDICTOR_HORIZONTAL . IMWRITE_JPEG2000_COMPRESSION_X1000 = 272,//!< For JPEG2000, use to specify the target compression rate (multiplied by 1000). The value can be from 0 to 1000. Default is 1000. @@ -540,8 +540,11 @@ can be saved using this function, with these exceptions: 32-bit signed (CV_32S), 32-bit float (CV_32F) and 64-bit float (CV_64F) images can be saved. - Multiple images (vector of Mat) can be saved in TIFF format (see the code sample below). - - 32-bit float 3-channel (CV_32FC3) TIFF images will be saved - using the LogLuv high dynamic range encoding (4 bytes per pixel) + - 32-bit float 3-channel (CV_32FC3) TIFF images can be saved + using the LogLuv high dynamic range encoding (4 bytes per pixel) through TIFF_COMPRESSION_SGILOG or + (3 bytes per pixel) through TIFF_COMPRESSION_SGILOG24. + - Other compression schemes (LZW...) are supported as well for 32F depth, but the efficiency might not + be very good for the floating-point representation bit patterns. - With GIF encoder, 8-bit unsigned (CV_8U) images can be saved. - GIF images with an alpha channel can be saved using this function. To achieve this, create an 8-bit 4-channel (CV_8UC4) BGRA image, ensuring the alpha channel is the last component. diff --git a/modules/imgcodecs/src/grfmt_tiff.cpp b/modules/imgcodecs/src/grfmt_tiff.cpp index 902e061e31..f47e4722bb 100644 --- a/modules/imgcodecs/src/grfmt_tiff.cpp +++ b/modules/imgcodecs/src/grfmt_tiff.cpp @@ -1312,8 +1312,12 @@ bool TiffEncoder::writeLibTiff( const std::vector& img_vec, const std::vect cv::Ptr tif_cleanup(tif, cv_tiffCloseHandle); //Settings that matter to all images - int compression = IMWRITE_TIFF_COMPRESSION_LZW; - int predictor = IMWRITE_TIFF_PREDICTOR_HORIZONTAL; + const int compression_default_32F = IMWRITE_TIFF_COMPRESSION_NONE; + const int compression_default = IMWRITE_TIFF_COMPRESSION_LZW; + const int predictor_default_32F = IMWRITE_TIFF_PREDICTOR_FLOATINGPOINT; + const int predictor_default = IMWRITE_TIFF_PREDICTOR_HORIZONTAL; + int compression = -1; + int predictor = -1; int resUnit = -1, dpiX = -1, dpiY = -1; if(readParam(params, IMWRITE_TIFF_COMPRESSION, compression)) @@ -1420,14 +1424,37 @@ bool TiffEncoder::writeLibTiff( const std::vector& img_vec, const std::vect CV_TIFF_CHECK_CALL(TIFFSetField(tif, TIFFTAG_PAGENUMBER, page, img_vec.size())); } - if (type == CV_32FC3 && compression == COMPRESSION_SGILOG) - { - if (!write_32FC3_SGILOG(img, tif)) - return false; - continue; - } + const bool is32F = (depth == CV_32F); + int page_compression = + (compression < 0) ? + (is32F ? compression_default_32F : compression_default) : + compression; + int page_predictor = + (predictor < 0) ? + (is32F ? predictor_default_32F : predictor_default) : + predictor; - int page_compression = compression; + if ((page_compression == COMPRESSION_SGILOG) || (page_compression == COMPRESSION_SGILOG24)) + { + if (depth != CV_32F) + CV_Error(cv::Error::StsError, "SGILOG requires 32F"); + else if ((page_compression == COMPRESSION_SGILOG24) && (type != CV_32FC3)) + CV_Error(cv::Error::StsError, "SGILOG24 requires 32FC3"); + else if ((page_compression == COMPRESSION_SGILOG) && (type != CV_32FC1) && (type != CV_32FC3)) + CV_Error(cv::Error::StsError, "SGILOG requires 32FC1 or 32FC3"); + else if ((page_compression == COMPRESSION_SGILOG) && (type == CV_32FC3)) + { + if (!write_32FC3_SGILOG(img, tif)) + return false; + continue; + } + else + { + if (!write_32F_SGILOG(img, tif, page_compression)) + return false; + continue; + } + } int bitsPerChannel = -1; uint16_t sample_format = SAMPLEFORMAT_INT; @@ -1460,7 +1487,6 @@ bool TiffEncoder::writeLibTiff( const std::vector& img_vec, const std::vect case CV_32F: { bitsPerChannel = 32; - page_compression = COMPRESSION_NONE; sample_format = SAMPLEFORMAT_IEEEFP; break; } @@ -1498,7 +1524,7 @@ bool TiffEncoder::writeLibTiff( const std::vector& img_vec, const std::vect if (page_compression == COMPRESSION_LZW || page_compression == COMPRESSION_ADOBE_DEFLATE || page_compression == COMPRESSION_DEFLATE) { - CV_TIFF_CHECK_CALL(TIFFSetField(tif, TIFFTAG_PREDICTOR, predictor)); + CV_TIFF_CHECK_CALL(TIFFSetField(tif, TIFFTAG_PREDICTOR, page_predictor)); } if (resUnit >= RESUNIT_NONE && resUnit <= RESUNIT_CENTIMETER) @@ -1583,6 +1609,41 @@ bool TiffEncoder::write_32FC3_SGILOG(const Mat& _img, void* tif_) return true; } +bool TiffEncoder::write_32F_SGILOG(const Mat& _img, void* tif_, int compression) +{ + TIFF* tif = (TIFF*)tif_; + CV_Assert(tif); + const int nChannels = _img.channels(); + CV_Assert( + ((compression == COMPRESSION_SGILOG) && ((nChannels == 1) || (nChannels == 3))) || + ((compression == COMPRESSION_SGILOG24) && (nChannels == 3)) + ); + + Mat img; + if (nChannels == 1) + img = _img; + else + cvtColor(_img, img, COLOR_BGR2XYZ); + + //done by caller: CV_TIFF_CHECK_CALL(TIFFSetField(tif, TIFFTAG_IMAGEWIDTH, img.cols)); + //done by caller: CV_TIFF_CHECK_CALL(TIFFSetField(tif, TIFFTAG_IMAGELENGTH, img.rows)); + CV_TIFF_CHECK_CALL(TIFFSetField(tif, TIFFTAG_SAMPLESPERPIXEL, nChannels)); + CV_TIFF_CHECK_CALL(TIFFSetField(tif, TIFFTAG_BITSPERSAMPLE, 32)); + CV_TIFF_CHECK_CALL(TIFFSetField(tif, TIFFTAG_COMPRESSION, compression)); + CV_TIFF_CHECK_CALL(TIFFSetField(tif, TIFFTAG_PHOTOMETRIC, + (nChannels == 1) ? PHOTOMETRIC_LOGL : PHOTOMETRIC_LOGLUV)); + CV_TIFF_CHECK_CALL(TIFFSetField(tif, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG)); + CV_TIFF_CHECK_CALL(TIFFSetField(tif, TIFFTAG_SGILOGDATAFMT, SGILOGDATAFMT_FLOAT)); + CV_TIFF_CHECK_CALL(TIFFSetField(tif, TIFFTAG_ROWSPERSTRIP, 1)); + const int strip_size = nChannels * img.cols; + for (int i = 0; i < img.rows; i++) + { + CV_TIFF_CHECK_CALL(TIFFWriteEncodedStrip(tif, i, (tdata_t)img.ptr(i), strip_size * sizeof(float)) != (tsize_t)-1); + } + CV_TIFF_CHECK_CALL(TIFFWriteDirectory(tif)); + return true; +} + bool TiffEncoder::writemulti(const std::vector& img_vec, const std::vector& params) { return writeLibTiff(img_vec, params); diff --git a/modules/imgcodecs/src/grfmt_tiff.hpp b/modules/imgcodecs/src/grfmt_tiff.hpp index 0d1f511372..d1dd0e3392 100644 --- a/modules/imgcodecs/src/grfmt_tiff.hpp +++ b/modules/imgcodecs/src/grfmt_tiff.hpp @@ -134,6 +134,7 @@ public: protected: bool writeLibTiff( const std::vector& img_vec, const std::vector& params ); bool write_32FC3_SGILOG(const Mat& img, void* tif); + bool write_32F_SGILOG(const Mat& img, void* tif, int compression); private: TiffEncoder(const TiffEncoder &); // copy disabled diff --git a/modules/imgcodecs/test/test_tiff.cpp b/modules/imgcodecs/test/test_tiff.cpp index af995158b4..400e9a9a15 100644 --- a/modules/imgcodecs/test/test_tiff.cpp +++ b/modules/imgcodecs/test/test_tiff.cpp @@ -934,7 +934,7 @@ Imgcodes_Tiff_TypeAndComp all_types[] = { { CV_16UC1, true }, { CV_16UC3, true }, { CV_16UC4, true }, { CV_16SC1, true }, { CV_16SC3, true }, { CV_16SC4, true }, { CV_32SC1, true }, { CV_32SC3, true }, { CV_32SC4, true }, - { CV_32FC1, false }, { CV_32FC3, false }, { CV_32FC4, false }, // No compression + { CV_32FC1, true }, { CV_32FC3, true }, { CV_32FC4, true }, { CV_64FC1, false }, { CV_64FC3, false }, { CV_64FC4, false } // No compression }; @@ -1293,6 +1293,47 @@ TEST(Imgcodecs_Tiff, read_junk) { ASSERT_TRUE(img.empty()); } + +typedef int Imgcodecs_Tiff_32F_Compressions_32F_Values; +typedef testing::TestWithParam Imgcodecs_Tiff_32F_Compressions_32F; + +TEST_P(Imgcodecs_Tiff_32F_Compressions_32F, compressions_32F) +{ + const int compression = GetParam(); + + const Size size(64, 64); + Mat src = Mat(size, CV_32FC1); + cv::randu(src, cv::Scalar::all(0.), cv::Scalar::all(1.)); + + std::vector params; + if (compression > 0) + { + params.push_back(IMWRITE_TIFF_COMPRESSION); + params.push_back(compression); + } + + std::vector encoded_data; + imencode(".tiff", src, encoded_data, params); + + Mat dst; + imdecode(encoded_data, IMREAD_UNCHANGED, &dst); + + EXPECT_LE(cvtest::norm(src, dst, NORM_INF), 1e-6); +} + +const int Imgcodecs_Tiff_32F_Compressions_32F_All_Values[] = +{ + -1,//will mean "default" + IMWRITE_TIFF_COMPRESSION_NONE, + IMWRITE_TIFF_COMPRESSION_LZW, + //IMWRITE_TIFF_COMPRESSION_LZMA,//might not be configured + //IMWRITE_TIFF_COMPRESSION_ZSTD,//might not be configured + //IMWRITE_TIFF_COMPRESSION_DEFLATE,//deprecated + IMWRITE_TIFF_COMPRESSION_ADOBE_DEFLATE, +}; + +INSTANTIATE_TEST_CASE_P(compressions_32F, Imgcodecs_Tiff_32F_Compressions_32F, testing::ValuesIn(Imgcodecs_Tiff_32F_Compressions_32F_All_Values)); + #endif }} // namespace From ae6198e0003ced466af6c6fdbb137b7fb8d219a8 Mon Sep 17 00:00:00 2001 From: Vincent Rabaud Date: Wed, 15 Apr 2026 12:43:47 +0200 Subject: [PATCH 08/19] Rewrite getBitsFromByteList to avoid harmless buffer overflow At the end, currentByte = byteList.ptr()[base + currentByteIdx]; could be out of bound though never used. --- .../objdetect/src/aruco/aruco_dictionary.cpp | 42 +++++++++---------- 1 file changed, 19 insertions(+), 23 deletions(-) diff --git a/modules/objdetect/src/aruco/aruco_dictionary.cpp b/modules/objdetect/src/aruco/aruco_dictionary.cpp index 7dfdcc3009..abacd9e4d1 100644 --- a/modules/objdetect/src/aruco/aruco_dictionary.cpp +++ b/modules/objdetect/src/aruco/aruco_dictionary.cpp @@ -235,34 +235,30 @@ Mat Dictionary::getBitsFromByteList(const Mat &byteList, int markerSize, int rot CV_Assert(byteList.channels() >= 4); CV_Assert(rotationId >= 0 && rotationId < 4); - Mat bits(markerSize, markerSize, CV_8UC1, Scalar::all(0)); - - unsigned char base2List[] = { 128, 64, 32, 16, 8, 4, 2, 1 }; + Mat bits = Mat::zeros(markerSize, markerSize, CV_8UC1); + unsigned char *bitsPtr = bits.ptr(); // Use a base offset for the selected rotation int nbytes = (bits.cols * bits.rows + 8 - 1) / 8; // integer ceil int base = rotationId * nbytes; - int currentByteIdx = 0; - unsigned char currentByte = byteList.ptr()[base + currentByteIdx]; - int currentBit = 0; + const unsigned char *currentBytePtr = byteList.ptr() + base; + const unsigned char *currentBytePtrEnd = currentBytePtr + bits.total() / 8; - for(int row = 0; row < bits.rows; row++) { - for(int col = 0; col < bits.cols; col++) { - if(currentByte >= base2List[currentBit]) { - bits.at(row, col) = 1; - currentByte -= base2List[currentBit]; - } - currentBit++; - if(currentBit == 8) { - currentByteIdx++; - currentByte = byteList.ptr()[base + currentByteIdx]; - // if not enough bits for one more byte, we are in the end - // update bit position accordingly - if(8 * (currentByteIdx + 1) > (int)bits.total()) - currentBit = 8 * (currentByteIdx + 1) - (int)bits.total(); - else - currentBit = 0; // ok, bits enough for next byte - } + for(;currentBytePtr < currentBytePtrEnd; ++currentBytePtr) { + unsigned char currentByte = *currentBytePtr; + for(int mask = 1 << 7; mask != 0; mask >>= 1) { + if (currentByte & mask) *bitsPtr = 1; + ++bitsPtr; + } + } + // if not enough bits for one more byte, we are in the end + // update bit position accordingly + if (bits.total() % 8 != 0) { + unsigned char currentByte = *currentBytePtrEnd; + int mask = 1 << ((bits.total() % 8) - 1); + for(; mask != 0; mask >>= 1) { + if (currentByte & mask) *bitsPtr = 1; + ++bitsPtr; } } return bits; From 3215d7e6ea9681d8e735eb77aa807f91d628102c Mon Sep 17 00:00:00 2001 From: vrooomy Date: Thu, 16 Apr 2026 18:09:36 +0530 Subject: [PATCH 09/19] fix Flatten axis=rank bug --- modules/dnn/src/onnx/onnx_importer.cpp | 40 +++++++++++++++++++++---- modules/dnn/test/test_onnx_importer.cpp | 1 + 2 files changed, 36 insertions(+), 5 deletions(-) diff --git a/modules/dnn/src/onnx/onnx_importer.cpp b/modules/dnn/src/onnx/onnx_importer.cpp index 11d24eadb8..11fde4278a 100644 --- a/modules/dnn/src/onnx/onnx_importer.cpp +++ b/modules/dnn/src/onnx/onnx_importer.cpp @@ -2225,7 +2225,9 @@ void ONNXImporter::parseFlatten(LayerParams& layerParams, const opencv_onnx::Nod { constBlobsExtraInfo.insert(std::make_pair(node_proto.output(0), getBlobExtraInfo(node_proto, 0))); } - int axis = normalize_axis(axis_, input.dims); + int axis = axis_; + if (axis < 0) axis += input.dims; + axis = std::max(0, std::min(axis, input.dims)); int out_size[2] = {1, 1}; for (int i = 0; i < axis; ++i) @@ -2244,18 +2246,46 @@ void ONNXImporter::parseFlatten(LayerParams& layerParams, const opencv_onnx::Nod IterShape_t shapeIt = outShapes.find(node_proto.input(0)); CV_Assert(shapeIt != outShapes.end()); MatShape inpShape = shapeIt->second; - int axis = normalize_axis(axis_, inpShape.size()); + int axis = axis_; + if (axis < 0) axis += (int)inpShape.size(); + axis = std::max(0, std::min(axis, (int)inpShape.size())); - if (axis == 0 || axis == inpShape.size()) + if (axis == (int)inpShape.size()) { LayerParams reshapeLp; reshapeLp.name = layerParams.name + "/reshape"; reshapeLp.type = "Reshape"; CV_Assert(layer_id.find(reshapeLp.name) == layer_id.end()); - - inpShape.insert(axis == 0 ? inpShape.begin() : inpShape.end(), 1); + inpShape.push_back(1); reshapeLp.set("dim", DictValue::arrayInt(&inpShape[0], inpShape.size())); + opencv_onnx::NodeProto proto; + proto.add_input(node_proto.input(0)); + proto.add_output(reshapeLp.name); + addLayer(reshapeLp, proto); + LayerParams flatLp; + flatLp.name = layerParams.name + "/flatten"; + flatLp.type = "Flatten"; + CV_Assert(layer_id.find(flatLp.name) == layer_id.end()); + flatLp.set("axis", 0); + flatLp.set("end_axis", (int)inpShape.size() - 2); + opencv_onnx::NodeProto proto2; + proto2.add_input(reshapeLp.name); + proto2.add_output(flatLp.name); + addLayer(flatLp, proto2); + layerParams.type = "Identity"; + node_proto.set_input(0, flatLp.name); + addLayer(layerParams, node_proto); + return; + } + if (axis == 0) + { + LayerParams reshapeLp; + reshapeLp.name = layerParams.name + "/reshape"; + reshapeLp.type = "Reshape"; + CV_Assert(layer_id.find(reshapeLp.name) == layer_id.end()); + inpShape.insert(inpShape.begin(), 1); + reshapeLp.set("dim", DictValue::arrayInt(&inpShape[0], inpShape.size())); opencv_onnx::NodeProto proto; proto.add_input(node_proto.input(0)); proto.add_output(reshapeLp.name); diff --git a/modules/dnn/test/test_onnx_importer.cpp b/modules/dnn/test/test_onnx_importer.cpp index a69e94990b..cd1ae9d2fd 100644 --- a/modules/dnn/test/test_onnx_importer.cpp +++ b/modules/dnn/test/test_onnx_importer.cpp @@ -1170,6 +1170,7 @@ TEST_P(Test_ONNX_layers, DynamicReshape) testONNXModels("dynamic_reshape_opset_11"); testONNXModels("flatten_by_prod"); testONNXModels("flatten_const"); + testONNXModels("flatten_axis_numaxes"); } TEST_P(Test_ONNX_layers, Reshape) From 166c235603f88acc84d2934269841a25e00c270a Mon Sep 17 00:00:00 2001 From: Jesus Armando Anaya Date: Fri, 17 Apr 2026 00:53:45 -0700 Subject: [PATCH 10/19] Merge pull request #28745 from JesusAnaya:fix-bundled-protobuf-include-priority build: fix bundled protobuf headers being overridden by system-installed protobuf #28745 On macOS with Apple Clang 17, the compiler implicitly injects `/usr/local/include` as a high-priority user include (`-I`). The `SYSTEM` keyword in `target_include_directories` emits `-isystem` for the bundled protobuf path, which is searched after user `-I` paths, allowing a system-installed protobuf v4+ (which requires C++17 and abseil-cpp) to override the bundled headers. Removing `SYSTEM` causes CMake to emit `-I` instead of `-isystem`, placing the bundled path before the implicit system include, restoring the intended header resolution order regardless of what protobuf version is installed on the host. ### Background Apple Clang 17 (shipped with macOS 26 Tahoe Command Line Tools) changed the compiler driver behavior: it now unconditionally adds `-I/usr/local/include` to the search path, which is where Homebrew installs headers. Previous Apple Clang versions either omitted this path or added it as a lower-priority system include. When Homebrew's `protobuf` is installed (v4+, e.g. `33.4_1`), its headers at `/usr/local/include/google/protobuf/` are resolved before the bundled `3rdparty/protobuf/src/` because `SYSTEM PUBLIC` in `target_include_directories` causes the bundled path to be emitted as `-isystem`, which ranks below `-I` in compiler search order. This was confirmed by inspecting the generated `flags.make`: ``` # before fix -isystem /path/to/opencv/3rdparty/protobuf/src # after fix -I /path/to/opencv/3rdparty/protobuf/src (first in the include list) ``` The system protobuf v4+ headers require C++17 and abseil-cpp. Since OpenCV 4.x compiles with `-std=c++11`, every translation unit in `3rdparty/protobuf/` fails: ``` /usr/local/include/google/protobuf/port_def.inc: error: Protobuf only supports C++17 and newer. /usr/local/include/absl/base/policy_checks.h: error: C++ versions less than C++17 are not supported. ``` Note that `include_directories(BEFORE ...)` on the line above was intended to prevent exactly this, but it was overridden by `target_include_directories(SYSTEM PUBLIC ...)` deduplicating the path into a single `-isystem` entry. This is essentially the same class of problem tracked in #13328. ### Impact The change is a one-line diff. Removing `SYSTEM` has no effect on Linux or Windows since those toolchains do not inject `/usr/local/include` implicitly. Warning suppression for protobuf headers in dependent modules is unaffected because OpenCV already applies explicit `-Wno-*` flags for all protobuf-related warnings in `3rdparty/protobuf/CMakeLists.txt`. ### Tested on - macOS 26.3 (Tahoe), Apple Clang 17.0.0 (`clang-1700.6.4.2`), Homebrew `protobuf 33.4_1`, full build passes with `BUILD_PROTOBUF=ON`. Fixes #28743, related to #27886. ### 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 - [x] There is a reference to the original bug report and related work - [ ] There is accuracy test, performance test and test data in opencv_extra repository, if applicable Patch to opencv_extra has the same branch name. - [ ] The feature is well documented and sample code can be built with the project CMake --- cmake/OpenCVCompilerOptions.cmake | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/cmake/OpenCVCompilerOptions.cmake b/cmake/OpenCVCompilerOptions.cmake index 512e6d1abd..5a0cb3d1bc 100644 --- a/cmake/OpenCVCompilerOptions.cmake +++ b/cmake/OpenCVCompilerOptions.cmake @@ -477,6 +477,14 @@ if(APPLE AND NOT CMAKE_CROSSCOMPILING AND NOT DEFINED ENV{LDFLAGS} AND EXISTS "/ link_directories("/usr/local/lib") endif() +if(APPLE AND NOT CMAKE_CROSSCOMPILING AND CV_CLANG AND EXISTS "/usr/local/include") + # Apple Clang 17+ implicitly injects -I/usr/local/include as a high-priority + # user include, causing system-installed headers (e.g. Homebrew protobuf v4+) + # to override bundled third-party libraries added via -isystem. + # Demote /usr/local/include to -isystem so bundled -isystem paths are searched first. + add_compile_options("-isystem/usr/local/include") +endif() + if(ENABLE_BUILD_HARDENING) include("${CMAKE_CURRENT_LIST_DIR}/OpenCVCompilerDefenses.cmake") endif() From 6aba6fda2d6f76250b599c46d11cc3b25e3001cc Mon Sep 17 00:00:00 2001 From: Alexander Smorkalov Date: Fri, 17 Apr 2026 10:56:46 +0300 Subject: [PATCH 11/19] Disable FFmpeg tests that need FFmpeg wrapper rebuild on Windows. --- modules/videoio/test/test_ffmpeg.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/modules/videoio/test/test_ffmpeg.cpp b/modules/videoio/test/test_ffmpeg.cpp index 0854dec433..986dbdb50a 100644 --- a/modules/videoio/test/test_ffmpeg.cpp +++ b/modules/videoio/test/test_ffmpeg.cpp @@ -845,6 +845,8 @@ TEST(videoio_ffmpeg, create_with_property_badarg) EXPECT_FALSE(cap.isOpened()); } +// requires FFmpeg wrapper rebuild on Windows +#ifndef _WIN32 TEST(videoio_ffmpeg, open_with_format_cv8uc3) { if (!videoio_registry::hasBackend(CAP_FFMPEG)) @@ -860,6 +862,7 @@ TEST(videoio_ffmpeg, open_with_format_cv8uc3) ASSERT_TRUE(cap.read(frame)); EXPECT_EQ(frame.channels(), 3); } +#endif // related issue: https://github.com/opencv/opencv/issues/16821 TEST(videoio_ffmpeg, DISABLED_open_from_web) From ab0def4854de69392b92afd3956cf7d877dce9b3 Mon Sep 17 00:00:00 2001 From: Anand Mahesh Date: Sat, 18 Apr 2026 13:29:41 +0530 Subject: [PATCH 12/19] chore/faster-ocl-bfmatcher-crosscheck The existing 2-pass cross-check implementation requires both forward and reverse distance computations on GPU, followed by CPU-side filtering. This adds a single-pass kernel using int64 atomics to eliminate the reverse pass. The BruteForceMatch_CrossCheckMatch kernel performs forward matching and tracks best train->query mappings in one GPU pass, reducing memory transfers. Falls back to 2-pass approach on devices without cl_khr_int64_base_atomics extension. Tested on NVIDIA RTX 5060 Ti: 1.65x-1.97x faster than 2-pass GPU, 3.49x-7.33x faster than CPU at 640x480 to 1920x1080 resolutions. --- modules/features2d/src/matchers.cpp | 125 ++++++++++++++-- .../src/opencl/brute_force_match.cl | 136 +++++++++++++++++- 2 files changed, 248 insertions(+), 13 deletions(-) diff --git a/modules/features2d/src/matchers.cpp b/modules/features2d/src/matchers.cpp index 950fd96d58..71026dbb15 100644 --- a/modules/features2d/src/matchers.cpp +++ b/modules/features2d/src/matchers.cpp @@ -74,7 +74,7 @@ static void ensureSizeIsEnough(int rows, int cols, int type, UMat &m) } static bool ocl_matchSingle(InputArray query, InputArray train, - UMat &trainIdx, UMat &distance, int distType) + UMat &trainIdx, UMat *distance, int distType) { if (query.empty() || train.empty()) return false; @@ -83,7 +83,8 @@ static bool ocl_matchSingle(InputArray query, InputArray train, const int query_cols = query.cols(); ensureSizeIsEnough(1, query_rows, CV_32S, trainIdx); - ensureSizeIsEnough(1, query_rows, CV_32F, distance); + if (distance) + ensureSizeIsEnough(1, query_rows, CV_32F, *distance); ocl::Device devDef = ocl::Device::getDefault(); @@ -117,7 +118,10 @@ static bool ocl_matchSingle(InputArray query, InputArray train, idx = k.set(idx, ocl::KernelArg::PtrReadOnly(uquery)); idx = k.set(idx, ocl::KernelArg::PtrReadOnly(utrain)); idx = k.set(idx, ocl::KernelArg::PtrWriteOnly(trainIdx)); - idx = k.set(idx, ocl::KernelArg::PtrWriteOnly(distance)); + UMat dummyDist; + if (!distance) + ensureSizeIsEnough(1, query_rows, CV_32F, dummyDist); + idx = k.set(idx, ocl::KernelArg::PtrWriteOnly(distance ? *distance : dummyDist)); idx = k.set(idx, uquery.rows); idx = k.set(idx, uquery.cols); idx = k.set(idx, utrain.rows); @@ -734,7 +738,7 @@ Ptr BFMatcher::clone( bool emptyTrainData ) const static bool ocl_match(InputArray query, InputArray _train, std::vector< std::vector > &matches, int dstType) { UMat trainIdx, distance; - if (!ocl_matchSingle(query, _train, trainIdx, distance, dstType)) + if (!ocl_matchSingle(query, _train, trainIdx, &distance, dstType)) return false; if (!ocl_matchDownload(trainIdx, distance, matches)) return false; @@ -753,22 +757,119 @@ static bool ocl_knnMatch(InputArray query, InputArray _train, std::vector< std:: return true; } -// Run match in both directions on GPU, cross-check filter on CPU. -// Reuses ocl_matchSingle so no new kernel is needed. +static bool ocl_matchWithCrossCheckSinglePass(InputArray query, InputArray train, + std::vector< std::vector >& matches, int dstType) +{ + if (query.empty() || train.empty()) + return false; + + const int query_rows = query.rows(); + const int train_rows = train.rows(); + + UMat fwdIdx, fwdDist; + ensureSizeIsEnough(1, query_rows, CV_32S, fwdIdx); + ensureSizeIsEnough(1, query_rows, CV_32F, fwdDist); + + UMat revBest; + ensureSizeIsEnough(1, train_rows, CV_32SC2, revBest); + revBest.setTo(Scalar::all(-1)); + + ocl::Device devDef = ocl::Device::getDefault(); + + UMat uquery = query.getUMat(), utrain = train.getUMat(); + int kercn = 1; + if (devDef.isIntel() && + (0 == (uquery.step % 4)) && (0 == (uquery.cols % 4)) && (0 == (uquery.offset % 4)) && + (0 == (utrain.step % 4)) && (0 == (utrain.cols % 4)) && (0 == (utrain.offset % 4))) + kercn = 4; + + int block_size = 16; + int max_desc_len = 0; + bool is_cpu = devDef.type() == ocl::Device::TYPE_CPU; + if (uquery.cols <= 64) + max_desc_len = 64 / kercn; + else if (uquery.cols <= 128 && !is_cpu) + max_desc_len = 128 / kercn; + + int depth = uquery.depth(); + cv::String opts; + opts = cv::format("-D T=%s -D TN=%s -D kercn=%d %s -D DIST_TYPE=%d -D BLOCK_SIZE=%d -D MAX_DESC_LEN=%d -D HAVE_INT64_ATOMICS", + ocl::typeToStr(depth), ocl::typeToStr(CV_MAKETYPE(depth, kercn)), kercn, depth == CV_32F ? "-D T_FLOAT" : "", dstType, block_size, max_desc_len); + ocl::Kernel k("BruteForceMatch_CrossCheckMatch", ocl::features2d::brute_force_match_oclsrc, opts); + if(k.empty()) + return false; + + size_t globalSize[] = {((size_t)uquery.size().height + block_size - 1) / block_size * block_size, (size_t)block_size}; + size_t localSize[] = {(size_t)block_size, (size_t)block_size}; + + int idx = 0; + idx = k.set(idx, ocl::KernelArg::PtrReadOnly(uquery)); + idx = k.set(idx, ocl::KernelArg::PtrReadOnly(utrain)); + idx = k.set(idx, ocl::KernelArg::PtrWriteOnly(fwdIdx)); + idx = k.set(idx, ocl::KernelArg::PtrWriteOnly(fwdDist)); + idx = k.set(idx, ocl::KernelArg::PtrWriteOnly(revBest)); + idx = k.set(idx, uquery.rows); + idx = k.set(idx, uquery.cols); + idx = k.set(idx, utrain.rows); + idx = k.set(idx, utrain.cols); + idx = k.set(idx, (int)(uquery.step / sizeof(float))); + + if (!k.run(2, globalSize, localSize, false)) + return false; + + Mat fwdIdxCPU = fwdIdx.getMat(ACCESS_READ); + Mat fwdDistCPU = fwdDist.getMat(ACCESS_READ); + Mat revBestCPU = revBest.getMat(ACCESS_READ); + + if (fwdIdxCPU.empty() || fwdDistCPU.empty() || revBestCPU.empty()) + return false; + + const int* fwd = fwdIdxCPU.ptr(); + const float* dist = fwdDistCPU.ptr(); + const int* revBestData = revBestCPU.ptr(); + + matches.clear(); + matches.reserve(query_rows); + + for (int q = 0; q < query_rows; ++q) + { + int t = fwd[q]; + if (t >= 0 && t < train_rows) + { + uint64_t packed_val; + memcpy(&packed_val, &revBestData[2 * t], sizeof(uint64_t)); + int revQueryIdx = (int)(packed_val & 0xFFFFFFFF); + if (revQueryIdx == q) + { + matches.push_back(std::vector(1, DMatch(q, t, 0, dist[q]))); + } + } + } + return true; +} + static bool ocl_matchWithCrossCheck(InputArray query, InputArray train, std::vector< std::vector >& matches, int dstType) { - // Forward pass: for each query descriptor find nearest in train + if (query.empty() || train.empty()) + return false; + + ocl::Device devDef = ocl::Device::getDefault(); + + if (devDef.isExtensionSupported("cl_khr_int64_base_atomics")) + { + if (ocl_matchWithCrossCheckSinglePass(query, train, matches, dstType)) + return true; + } + UMat fwdIdx, fwdDist; - if (!ocl_matchSingle(query, train, fwdIdx, fwdDist, dstType)) + if (!ocl_matchSingle(query, train, fwdIdx, &fwdDist, dstType)) return false; - // Reverse pass: for each train descriptor find nearest in query - UMat revIdx, revDist; - if (!ocl_matchSingle(train, query, revIdx, revDist, dstType)) + UMat revIdx; + if (!ocl_matchSingle(train, query, revIdx, nullptr, dstType)) return false; - // Download index arrays (1 x N ints each — cheap) Mat fwdIdxCPU = fwdIdx.getMat(ACCESS_READ); Mat revIdxCPU = revIdx.getMat(ACCESS_READ); Mat fwdDistCPU = fwdDist.getMat(ACCESS_READ); diff --git a/modules/features2d/src/opencl/brute_force_match.cl b/modules/features2d/src/opencl/brute_force_match.cl index 8f0e183799..68f2590741 100644 --- a/modules/features2d/src/opencl/brute_force_match.cl +++ b/modules/features2d/src/opencl/brute_force_match.cl @@ -557,4 +557,138 @@ __kernel void BruteForceMatch_knnMatch( bestTrainIdx[queryIdx] = (int2)(myBestTrainIdx1, myBestTrainIdx2); bestDistance[queryIdx] = (float2)(myBestDistance1, myBestDistance2); } -} \ No newline at end of file +} + +#ifdef HAVE_INT64_ATOMICS +#pragma OPENCL EXTENSION cl_khr_int64_base_atomics : enable + +__kernel void BruteForceMatch_CrossCheckMatch( + __global T *query, + __global T *train, + __global int *bestTrainIdx, + __global float *bestDistance, + __global ulong *revBest, + int query_rows, + int query_cols, + int train_rows, + int train_cols, + int step +) +{ + const int lidx = get_local_id(0); + const int lidy = get_local_id(1); + const int groupidx = get_group_id(0); + + const int queryIdx = mad24(BLOCK_SIZE, groupidx, lidy); + const int queryOffset = min(queryIdx, query_rows - 1) * step; + __global TN *query_vec = (__global TN *)(query + queryOffset); + query_cols /= kercn; + + __local float sharebuffer[SHARED_MEM_SZ]; + __local value_type *s_query = (__local value_type *)sharebuffer; + +#if 0 < MAX_DESC_LEN + __local value_type *s_train = (__local value_type *)sharebuffer + BLOCK_SIZE * MAX_DESC_LEN; + #pragma unroll + for (int i = 0; i < MAX_DESC_LEN / BLOCK_SIZE; i++) + { + const int loadx = mad24(BLOCK_SIZE, i, lidx); + s_query[mad24(MAX_DESC_LEN, lidy, loadx)] = loadx < query_cols ? query_vec[loadx] : 0; + } +#else + __local value_type *s_train = (__local value_type *)sharebuffer + BLOCK_SIZE_ODD * BLOCK_SIZE; + const int s_query_i = mad24(BLOCK_SIZE_ODD, lidy, lidx); + const int s_train_i = mad24(BLOCK_SIZE_ODD, lidx, lidy); +#endif + + float myBestDistance = MAX_FLOAT; + int myBestTrainIdx = -1; + + for (int t = 0, endt = (train_rows + BLOCK_SIZE - 1) / BLOCK_SIZE; t < endt; t++) + { + result_type result = 0; + + const int trainOffset = min(mad24(BLOCK_SIZE, t, lidy), train_rows - 1) * step; + __global TN *train_vec = (__global TN *)(train + trainOffset); +#if 0 < MAX_DESC_LEN + #pragma unroll + for (int i = 0; i < MAX_DESC_LEN / BLOCK_SIZE; i++) + { + const int loadx = mad24(BLOCK_SIZE, i, lidx); + s_train[mad24(BLOCK_SIZE, lidx, lidy)] = loadx < train_cols ? train_vec[loadx] : 0; + + barrier(CLK_LOCAL_MEM_FENCE); + + result += reduce_multi_block(s_query, s_train, i, lidx, lidy); + + barrier(CLK_LOCAL_MEM_FENCE); + } +#else + for (int i = 0, endq = (query_cols + BLOCK_SIZE - 1) / BLOCK_SIZE; i < endq; i++) + { + const int loadx = mad24(BLOCK_SIZE, i, lidx); + if (loadx < query_cols) + { + s_query[s_query_i] = query_vec[loadx]; + s_train[s_train_i] = train_vec[loadx]; + } + else + { + s_query[s_query_i] = 0; + s_train[s_train_i] = 0; + } + + barrier(CLK_LOCAL_MEM_FENCE); + + result += reduce_block_match(s_query, s_train, lidx, lidy); + + barrier(CLK_LOCAL_MEM_FENCE); + } +#endif + result = DIST_RES(result); + + const int trainIdx = mad24(BLOCK_SIZE, t, lidx); + + if (queryIdx < query_rows && trainIdx < train_rows) + { + if (result < myBestDistance) + { + myBestDistance = result; + myBestTrainIdx = trainIdx; + } + + uint dist_bits = as_uint(result); + ulong packed = ((ulong)dist_bits << 32) | (ulong)queryIdx; + atom_min(&revBest[trainIdx], packed); + } + } + + barrier(CLK_LOCAL_MEM_FENCE); + + __local float *s_distance = (__local float *)sharebuffer; + __local int *s_trainIdx = (__local int *)(sharebuffer + BLOCK_SIZE_ODD * BLOCK_SIZE); + + s_distance += lidy * BLOCK_SIZE_ODD; + s_trainIdx += lidy * BLOCK_SIZE_ODD; + s_distance[lidx] = myBestDistance; + s_trainIdx[lidx] = myBestTrainIdx; + + barrier(CLK_LOCAL_MEM_FENCE); + + #pragma unroll + for (int k = 0; k < BLOCK_SIZE; k++) + { + if (myBestDistance > s_distance[k]) + { + myBestDistance = s_distance[k]; + myBestTrainIdx = s_trainIdx[k]; + } + } + + if (queryIdx < query_rows && lidx == 0) + { + bestTrainIdx[queryIdx] = myBestTrainIdx; + bestDistance[queryIdx] = myBestDistance; + } +} +#endif \ No newline at end of file From 1ae5a8586b3beb28864c79c6ea6ebc6e06ecfe94 Mon Sep 17 00:00:00 2001 From: Kumataro Date: Mon, 20 Apr 2026 02:11:33 +0900 Subject: [PATCH 13/19] Merge pull request #28615 from Kumataro:fix28606 imgcodecs(png): support cICP metadata for imreadWithMetadata() #28615 Close https://github.com/opencv/opencv/issues/28606 Related https://github.com/opencv/opencv/pull/27741 (support for imwriteWithMetadata) ### 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 - [x] 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. - [ ] The feature is well documented and sample code can be built with the project CMake --- CMakeLists.txt | 10 +++++++ modules/imgcodecs/CMakeLists.txt | 9 +++++++ modules/imgcodecs/src/grfmt_png.cpp | 14 ++++++++++ modules/imgcodecs/test/test_exif.cpp | 39 ++++++++++++++++++++++++++++ modules/imgcodecs/test/test_png.cpp | 11 +++++++- 5 files changed, 82 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 8e30518445..4138c81098 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1537,6 +1537,8 @@ if(WITH_SPNG) elseif(HAVE_SPNG) status(" PNG:" "${SPNG_LIBRARY} (ver ${SPNG_VERSION})") endif() + + status(" Metadata Support:" "EXIF XMP ICC") # SPNG does not support cICP chunk. elseif(WITH_PNG OR HAVE_PNG) status(" PNG:" PNG_FOUND THEN "${PNG_LIBRARY} (ver ${PNG_VERSION_STRING})" ELSE "build (ver ${PNG_VERSION_STRING})") if(BUILD_PNG AND PNG_HARDWARE_OPTIMIZATIONS) @@ -1565,6 +1567,14 @@ elseif(WITH_PNG OR HAVE_PNG) elseif(BUILD_PNG) status(" SIMD Support Request:" "NO") endif() + + if(NOT (PNG_VERSION_STRING VERSION_LESS "1.6.45")) + status(" Metadata Support:" "EXIF XMP ICC cICP") + elseif(NOT (PNG_VERSION_STRING VERSION_LESS "1.6.31")) + status(" Metadata Support:" "EXIF XMP ICC") + else() + status(" Metadata Support:" "XMP ICC") + endif() endif() if(WITH_TIFF OR HAVE_TIFF) diff --git a/modules/imgcodecs/CMakeLists.txt b/modules/imgcodecs/CMakeLists.txt index c322645b20..d11b26b9ca 100644 --- a/modules/imgcodecs/CMakeLists.txt +++ b/modules/imgcodecs/CMakeLists.txt @@ -204,4 +204,13 @@ if(TARGET opencv_test_imgcodecs AND ((HAVE_PNG AND NOT (PNG_VERSION_STRING VERSI # details: https://github.com/glennrp/libpng/commit/68cb0aaee3de6371b81a4613476d9b33e43e95b1 ocv_target_compile_definitions(opencv_test_imgcodecs PRIVATE OPENCV_IMGCODECS_PNG_WITH_EXIF=1) endif() +if(TARGET opencv_test_imgcodecs AND ((HAVE_PNG AND NOT (PNG_VERSION_STRING VERSION_LESS "1.6.45")) OR HAVE_SPNG)) + # cICP does not support in spng + # cICP support added in libpng 1.6.45 + ocv_target_compile_definitions(opencv_test_imgcodecs PRIVATE OPENCV_IMGCODECS_PNG_WITH_cICP=1) +endif() +if(TARGET opencv_test_imgcodecs AND HAVE_PNG AND PNG_VERSION_STRING VERSION_LESS "1.6.0") + # Old libpng (< 1.6.0) is known to have lower precision in internal RGB-to-Gray calculation for 16-bit images. + ocv_target_compile_definitions(opencv_test_imgcodecs PRIVATE OPENCV_IMGCODECS_PNG_EPS_16BIT_GRAY=9) +endif() ocv_add_perf_tests() diff --git a/modules/imgcodecs/src/grfmt_png.cpp b/modules/imgcodecs/src/grfmt_png.cpp index e1e5d9a710..824cf961f9 100644 --- a/modules/imgcodecs/src/grfmt_png.cpp +++ b/modules/imgcodecs/src/grfmt_png.cpp @@ -684,6 +684,20 @@ bool PngDecoder::readData( Mat& img ) m_exif.parseExif(exif, num_exif); } #endif +#ifdef PNG_cICP_SUPPORTED + png_byte prim_id, tran_id, matrix_id, video_full_range_flag; + if (png_get_cICP(m_png_ptr, m_info_ptr, &prim_id, &tran_id, &matrix_id, &video_full_range_flag)) + { + uint8_t cicp_data[4] = { + static_cast(prim_id), + static_cast(tran_id), + static_cast(matrix_id), + static_cast(video_full_range_flag) + }; + auto& out = m_metadata[IMAGE_METADATA_CICP]; + out.insert(out.end(), cicp_data, cicp_data + 4); + } +#endif result = true; } diff --git a/modules/imgcodecs/test/test_exif.cpp b/modules/imgcodecs/test/test_exif.cpp index bd89728ce0..ebc849e9e6 100644 --- a/modules/imgcodecs/test/test_exif.cpp +++ b/modules/imgcodecs/test/test_exif.cpp @@ -144,6 +144,17 @@ namespace opencv_test { namespace { return iccp_data; } +#ifdef OPENCV_IMGCODECS_PNG_WITH_cICP + static std::vector getSampleCicpData() { + return { + 9, // BT.2020 / BT.2100 + 16, // SMPTE ST 2084 (PQ) + 0, // Identity (RGB) + 1, // Full Range + }; + } +#endif + /** * Test to check whether the EXIF orientation tag was processed successfully or not. * The test uses a set of 8 images named testExifOrientation_{1 to 8}.(extension). @@ -457,18 +468,27 @@ TEST(Imgcodecs_Png, Read_Write_With_Exif) EXPECT_EQ(img2.rows, img.rows); EXPECT_EQ(img2.type(), imgtype); EXPECT_EQ(read_metadata_types, read_metadata_types2); + +#ifdef OPENCV_IMGCODECS_PNG_WITH_EXIF ASSERT_GE(read_metadata_types.size(), 1u); EXPECT_EQ(read_metadata, read_metadata2); EXPECT_EQ(read_metadata_types[0], IMAGE_METADATA_EXIF); EXPECT_EQ(read_metadata_types.size(), read_metadata.size()); EXPECT_EQ(read_metadata[0], metadata[0]); +#else + ASSERT_GE(read_metadata_types.size(), 0u); +#endif EXPECT_EQ(cv::norm(img2, img3, NORM_INF), 0.); double mse = cv::norm(img, img2, NORM_L2SQR)/(img.rows*img.cols); EXPECT_EQ(mse, 0); // png is lossless remove(outputname.c_str()); } +#ifdef OPENCV_IMGCODECS_PNG_WITH_cICP +TEST(Imgcodecs_Png, Read_Write_With_Exif_Xmp_Iccp_cICP) +#else TEST(Imgcodecs_Png, Read_Write_With_Exif_Xmp_Iccp) +#endif { int png_compression = 3; int imgtype = CV_MAKETYPE(CV_8U, 3); @@ -482,6 +502,11 @@ TEST(Imgcodecs_Png, Read_Write_With_Exif_Xmp_Iccp) getSampleIccpData(), }; +#ifdef OPENCV_IMGCODECS_PNG_WITH_cICP + metadata_types.push_back(IMAGE_METADATA_CICP); + metadata.push_back(getSampleCicpData()); +#endif + std::vector write_params = { IMWRITE_PNG_COMPRESSION, png_compression }; @@ -498,9 +523,23 @@ TEST(Imgcodecs_Png, Read_Write_With_Exif_Xmp_Iccp) EXPECT_EQ(img2.rows, img.rows); EXPECT_EQ(img2.type(), imgtype); +#ifdef OPENCV_IMGCODECS_PNG_WITH_EXIF EXPECT_EQ(metadata_types, read_metadata_types); EXPECT_EQ(read_metadata_types, read_metadata_types2); EXPECT_EQ(metadata, read_metadata); +#else + ASSERT_GE(read_metadata_types.size(), 2u); + EXPECT_EQ(read_metadata_types[0], IMAGE_METADATA_XMP); + EXPECT_EQ(read_metadata_types[1], IMAGE_METADATA_ICCP); + + ASSERT_GE(read_metadata_types2.size(), 2u); + EXPECT_EQ(read_metadata_types2[0], IMAGE_METADATA_XMP); + EXPECT_EQ(read_metadata_types2[1], IMAGE_METADATA_ICCP); + + ASSERT_GE(read_metadata.size(), 2u); + EXPECT_EQ(metadata[1], read_metadata[0]); + EXPECT_EQ(metadata[2], read_metadata[1]); +#endif remove(outputname.c_str()); } diff --git a/modules/imgcodecs/test/test_png.cpp b/modules/imgcodecs/test/test_png.cpp index 1c8eae233c..2110bcd795 100644 --- a/modules/imgcodecs/test/test_png.cpp +++ b/modules/imgcodecs/test/test_png.cpp @@ -8,6 +8,13 @@ namespace opencv_test { namespace { #if defined(HAVE_PNG) || defined(HAVE_SPNG) +// See https://github.com/opencv/opencv/pull/28615 +// Precision differences in 16-bit grayscale conversion between old and modern libpng versions +#define OPENCV_IMGCODECS_PNG_EPS_DEFAULT (4) +#ifndef OPENCV_IMGCODECS_PNG_EPS_16BIT_GRAY +#define OPENCV_IMGCODECS_PNG_EPS_16BIT_GRAY (OPENCV_IMGCODECS_PNG_EPS_DEFAULT) +#endif + TEST(Imgcodecs_Png, write_big) { const string root = cvtest::TS::ptr()->get_data_path(); @@ -276,10 +283,12 @@ TEST_P(Imgcodecs_Png_PngSuite, decode) cvtColor(gt_3, gt_258, COLOR_BGR2RGB); } + const double epsGrayAnydepth = ((gt.depth() == CV_16U) && (gt.channels() > 1)) ? OPENCV_IMGCODECS_PNG_EPS_16BIT_GRAY: OPENCV_IMGCODECS_PNG_EPS_DEFAULT; + // Perform comparisons with different imread flags EXPECT_PRED_FORMAT2(cvtest::MatComparator(1, 0), imread(filename, IMREAD_GRAYSCALE), gt_0); EXPECT_PRED_FORMAT2(cvtest::MatComparator(1, 0), imread(filename, IMREAD_COLOR), gt_1); - EXPECT_PRED_FORMAT2(cvtest::MatComparator(4, 0), imread(filename, IMREAD_ANYDEPTH), gt_2); + EXPECT_PRED_FORMAT2(cvtest::MatComparator(epsGrayAnydepth, 0), imread(filename, IMREAD_ANYDEPTH), gt_2); // IMREAD_GRAYSCALE is used. EXPECT_PRED_FORMAT2(cvtest::MatComparator(0, 0), imread(filename, IMREAD_COLOR | IMREAD_ANYDEPTH), gt_3); EXPECT_PRED_FORMAT2(cvtest::MatComparator(1, 0), imread(filename, IMREAD_COLOR_RGB), gt_256); EXPECT_PRED_FORMAT2(cvtest::MatComparator(0, 0), imread(filename, IMREAD_COLOR_RGB | IMREAD_ANYDEPTH), gt_258); From 30178be8008087d43946f97a46395b4f639c5012 Mon Sep 17 00:00:00 2001 From: Alex <117711477+0lekW@users.noreply.github.com> Date: Mon, 20 Apr 2026 23:58:58 +1200 Subject: [PATCH 14/19] Merge pull request #28833 from 0lekW:avfoundation-fractional-fps-seek-fix Avfoundation fractional fps seek fix #28833 Accompanied by PR on [opencv_extra](https://github.com/opencv/opencv_extra/pull/1345) Fixes https://github.com/opencv/opencv/issues/28831 AVFoundation backend seeking slowly drifts for any video which has a non-integer frame rate. Narrowing of a float `mAssetTrack.nominalFrameRate` to int32_t means fractional rates are treated as integer when seeking by frame number. Solution implemented in this branch is to use the tracks native rational frame duration when seeking (and available). I've also included an AVF only test case which uses a new test file in opencv_extra. This test just compares the video caps expected ms vs actual for a known rate file (23.976). Could consider extending this test to more then just AVF, as from what I can see there are no fractional rate tests / files. ### 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 - [x] 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/videoio/src/cap_avfoundation.mm | 16 +++++++-- modules/videoio/src/cap_avfoundation_mac.mm | 14 ++++++-- modules/videoio/test/test_video_io.cpp | 39 +++++++++++++++++++++ 3 files changed, 64 insertions(+), 5 deletions(-) diff --git a/modules/videoio/src/cap_avfoundation.mm b/modules/videoio/src/cap_avfoundation.mm index 4f74ddfaab..aa7ee57142 100644 --- a/modules/videoio/src/cap_avfoundation.mm +++ b/modules/videoio/src/cap_avfoundation.mm @@ -1104,9 +1104,19 @@ bool CvCaptureFile::setProperty(int property_id, double value) { t.value = value * t.timescale / 1000; retval = setupReadingAt(t); break; - case cv::CAP_PROP_POS_FRAMES: - retval = mAssetTrack.nominalFrameRate > 0 ? setupReadingAt(CMTimeMake(value, mAssetTrack.nominalFrameRate)) : false; + case cv::CAP_PROP_POS_FRAMES: { + // Build the seek CMTime from the track's native rational frame + // duration so non-integer rates (e.g. 23.976 = 24000/1001) are exact. + CMTime frameDuration = mAssetTrack.minFrameDuration; + if (frameDuration.timescale > 0 && frameDuration.value > 0) { + retval = setupReadingAt(CMTimeMultiply(frameDuration, (int32_t)value)); + } else if (mAssetTrack.nominalFrameRate > 0) { + retval = setupReadingAt(CMTimeMake(value, mAssetTrack.nominalFrameRate)); + } else { + retval = false; + } break; + } case cv::CAP_PROP_POS_AVI_RATIO: t = mAsset.duration; t.value = round(t.value * value); @@ -1349,4 +1359,4 @@ void CvVideoWriter_AVFoundation::write(cv::InputArray image) { } } -#pragma clang diagnostic pop \ No newline at end of file +#pragma clang diagnostic pop diff --git a/modules/videoio/src/cap_avfoundation_mac.mm b/modules/videoio/src/cap_avfoundation_mac.mm index 5ba1c01c60..3bcd0928de 100644 --- a/modules/videoio/src/cap_avfoundation_mac.mm +++ b/modules/videoio/src/cap_avfoundation_mac.mm @@ -1033,9 +1033,19 @@ bool CvCaptureFile::setProperty_(int property_id, double value) { t.value = value * t.timescale / 1000; retval = setupReadingAt(t); break; - case cv::CAP_PROP_POS_FRAMES: - retval = mAssetTrack.nominalFrameRate > 0 ? setupReadingAt(CMTimeMake(value, mAssetTrack.nominalFrameRate)) : false; + case cv::CAP_PROP_POS_FRAMES: { + // Build the seek CMTime from the track's native rational frame + // duration so non-integer rates (e.g. 23.976 = 24000/1001) are exact. + CMTime frameDuration = mAssetTrack.minFrameDuration; + if (frameDuration.timescale > 0 && frameDuration.value > 0) { + retval = setupReadingAt(CMTimeMultiply(frameDuration, (int32_t)value)); + } else if (mAssetTrack.nominalFrameRate > 0) { + retval = setupReadingAt(CMTimeMake(value, mAssetTrack.nominalFrameRate)); + } else { + retval = false; + } break; + } case cv::CAP_PROP_POS_AVI_RATIO: t = mAsset.duration; t.value = round(t.value * value); diff --git a/modules/videoio/test/test_video_io.cpp b/modules/videoio/test/test_video_io.cpp index 9c1b5cb209..784b89d600 100644 --- a/modules/videoio/test/test_video_io.cpp +++ b/modules/videoio/test/test_video_io.cpp @@ -1238,4 +1238,43 @@ inline static std::string stream_capture_ffmpeg_name_printer(const testing::Test INSTANTIATE_TEST_CASE_P(videoio, stream_capture_ffmpeg, testing::Values("h264", "h265", "mjpg.avi"), stream_capture_ffmpeg_name_printer); +// Seek by frame index must be accurate on non-integer frame rates. +typedef testing::TestWithParam PreciseSeekingTest; + +TEST_P(PreciseSeekingTest, seek_nonInteger_fps_frame_accurate) +{ + VideoCaptureAPIs apiPref = GetParam(); + if (!videoio_registry::hasBackend(apiPref)) + throw SkipTestException(cv::String("Backend is not available/disabled: ") + cv::videoio_registry::getBackendName(apiPref)); + + const std::string file = findDataFile("video/sample_23976fps.mp4"); + + VideoCapture cap; + ASSERT_TRUE(cap.open(file, apiPref)) + << "AVFoundation could not open " << file; + + const double fps_num = 24000.0; + const double fps_den = 1001.0; + const double frame_period_ms = fps_den * 1000.0 / fps_num; // ~41.7083 + const double tol_ms = frame_period_ms * 0.5; + + const int targets[] = { 0, 1, 24, 50, 99 }; + for (int N : targets) + { + ASSERT_TRUE(cap.set(CAP_PROP_POS_FRAMES, N)) + << "seek to frame " << N << " failed"; + Mat img; + ASSERT_TRUE(cap.read(img)) + << "read after seek to frame " << N << " failed"; + const double got_ms = cap.get(CAP_PROP_POS_MSEC); + const double expect_ms = N * frame_period_ms; + EXPECT_NEAR(got_ms, expect_ms, tol_ms) + << "non-integer fps seek drifted at frame N=" << N; + } +} + +VideoCaptureAPIs seekable_backeinds[] = {CAP_FFMPEG, CAP_MSMF, CAP_AVFOUNDATION}; + +INSTANTIATE_TEST_CASE_P(videoio, PreciseSeekingTest, testing::ValuesIn(seekable_backeinds), safe_capture_name_printer); + } // namespace From ba879cd60f505c3239a6bc5c4b8eb6ca87d143c2 Mon Sep 17 00:00:00 2001 From: jiajia Qian Date: Thu, 16 Apr 2026 10:14:05 +0800 Subject: [PATCH 15/19] core/ocl: fix incorrect results for in-place flip on strict OpenCL implementations MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Previously, OpenCL in-place flip kernels (rows/cols/both) could produce incorrect results compared to the CPU implementation on strict OpenCL drivers such as Mesa. The kernels relied on implicit load–load–store–store (LLSS) ordering when src and dst alias, which is not guaranteed by the OpenCL memory model and may be reordered. Some vendor drivers happened to preserve the expected ordering, masking the issue, but Mesa correctly exposes the undefined behavior. This change introduces dedicated in-place flip kernels that: - Explicitly detect in-place execution (src == dst) - Stage data through local memory tiles - Enforce correct ordering with work-group barriers - Avoid global memory read/write aliasing hazards The non in-place path is unchanged. With this fix, OpenCL in-place flip produces correct and consistent results across drivers, matches CPU behavior, and complies with the OpenCL memory model. Signed-off-by: jiajia Qian --- modules/core/src/matrix_transform.cpp | 71 +++++++--- modules/core/src/opencl/flip.cl | 196 ++++++++++++++++++++++++++ 2 files changed, 250 insertions(+), 17 deletions(-) diff --git a/modules/core/src/matrix_transform.cpp b/modules/core/src/matrix_transform.cpp index ce859e8eb8..245cd89843 100644 --- a/modules/core/src/matrix_transform.cpp +++ b/modules/core/src/matrix_transform.cpp @@ -973,7 +973,12 @@ static bool ocl_flip(InputArray _src, OutputArray _dst, int flipCode ) if (cn > 4) return false; - const char * kernelName; + Size size = _src.size(); + _dst.create(size, type); + UMat src = _src.getUMat(), dst = _dst.getUMat(); + bool inplace = (dst.u == src.u); + + String kernelName; if (flipCode == 0) kernelName = "arithm_flip_rows", flipType = FLIP_ROWS; else if (flipCode > 0) @@ -981,33 +986,65 @@ static bool ocl_flip(InputArray _src, OutputArray _dst, int flipCode ) else kernelName = "arithm_flip_rows_cols", flipType = FLIP_BOTH; + if(inplace) + kernelName += "_inplace"; + int pxPerWIy = (dev.isIntel() && (dev.type() & ocl::Device::TYPE_GPU)) ? 4 : 1; kercn = (cn!=3 || flipType == FLIP_ROWS) ? std::max(kercn, cn) : cn; + const int TILE_SIZE = 32, BLOCK_ROWS = 8; - ocl::Kernel k(kernelName, ocl::core::flip_oclsrc, - format( "-D T=%s -D T1=%s -D DEPTH=%d -D cn=%d -D PIX_PER_WI_Y=%d -D kercn=%d", + ocl::Kernel k(kernelName.c_str(), ocl::core::flip_oclsrc, + format( "-D T=%s -D T1=%s -D DEPTH=%d -D cn=%d -D PIX_PER_WI_Y=%d -D kercn=%d -D TILE_SIZE=%d -D BLOCK_ROWS=%d%s", kercn != cn ? ocl::typeToStr(CV_MAKE_TYPE(depth, kercn)) : ocl::vecopTypeToStr(CV_MAKE_TYPE(depth, kercn)), - kercn != cn ? ocl::typeToStr(depth) : ocl::vecopTypeToStr(depth), depth, cn, pxPerWIy, kercn)); + kercn != cn ? ocl::typeToStr(depth) : ocl::vecopTypeToStr(depth), depth, cn, pxPerWIy, kercn, TILE_SIZE, BLOCK_ROWS, + inplace ? " -D INPLACE" : "")); if (k.empty()) return false; - Size size = _src.size(); - _dst.create(size, type); - UMat src = _src.getUMat(), dst = _dst.getUMat(); - int cols = size.width * cn / kercn, rows = size.height; - cols = flipType == FLIP_COLS ? (cols + 1) >> 1 : cols; - rows = flipType & FLIP_ROWS ? (rows + 1) >> 1 : rows; + int work_cols = flipType == FLIP_COLS ? (cols + 1) >> 1 : cols; + int work_rows = flipType & FLIP_ROWS ? (rows + 1) >> 1 : rows; - k.args(ocl::KernelArg::ReadOnlyNoSize(src), - ocl::KernelArg::WriteOnly(dst, cn, kercn), rows, cols); + if (inplace) + { + k.args(ocl::KernelArg::ReadWriteNoSize(dst), rows, cols); - size_t maxWorkGroupSize = dev.maxWorkGroupSize(); - CV_Assert(maxWorkGroupSize % 4 == 0); + int gs_cols, gs_rows; + if (flipType == FLIP_COLS) + { + gs_cols = work_cols; + gs_rows = rows; + } + else if (flipType == FLIP_ROWS) + { + gs_cols = cols; + gs_rows = work_rows; + } + else // FLIP_BOTH + { + gs_cols = cols; + gs_rows = rows; + } - size_t globalsize[2] = { (size_t)cols, ((size_t)rows + pxPerWIy - 1) / pxPerWIy }, - localsize[2] = { maxWorkGroupSize / 4, 4 }; - return k.run(2, globalsize, (flipType == FLIP_COLS) && !dev.isIntel() ? localsize : NULL, false); + size_t globalsize[2] = { + (size_t)divUp(gs_cols, TILE_SIZE) * TILE_SIZE, + (size_t)divUp(gs_rows, TILE_SIZE) * BLOCK_ROWS + }; + size_t localsize[2] = { TILE_SIZE, BLOCK_ROWS }; + return k.run(2, globalsize, localsize, false); + } + else + { + k.args(ocl::KernelArg::ReadOnlyNoSize(src), + ocl::KernelArg::WriteOnly(dst, cn, kercn), work_rows, work_cols); + + size_t maxWorkGroupSize = dev.maxWorkGroupSize(); + CV_Assert(maxWorkGroupSize % 4 == 0); + + size_t globalsize[2] = { (size_t)work_cols, ((size_t)work_rows + pxPerWIy - 1) / pxPerWIy }; + size_t localsize[2] = { maxWorkGroupSize / 4, 4 }; + return k.run(2, globalsize, (flipType == FLIP_COLS) && !dev.isIntel() ? localsize : NULL, false); + } } #endif diff --git a/modules/core/src/opencl/flip.cl b/modules/core/src/opencl/flip.cl index afd14e4e1f..dbe23c73a4 100644 --- a/modules/core/src/opencl/flip.cl +++ b/modules/core/src/opencl/flip.cl @@ -63,6 +63,9 @@ #endif #define TSIZE ((int)sizeof(T1)*3) #endif +#define LDS_STEP (TILE_SIZE + 1) + +#ifndef INPLACE __kernel void arithm_flip_rows(__global const uchar * srcptr, int src_step, int src_offset, __global uchar * dstptr, int dst_step, int dst_offset, @@ -183,3 +186,196 @@ __kernel void arithm_flip_cols(__global const uchar * srcptr, int src_step, int } } } + +#else + +__kernel void arithm_flip_rows_inplace(__global uchar * srcptr, int src_step, int src_offset, + int rows, int cols) +{ + int gp_x = get_group_id(0); + int gp_y = get_group_id(1); + int lx = get_local_id(0); + int ly = get_local_id(1); + + __local T tile_top[TILE_SIZE * LDS_STEP]; + __local T tile_bottom[TILE_SIZE * LDS_STEP]; + + int half_rows = (rows + 1) / 2; + + int x = gp_x * TILE_SIZE + lx; + int y_top = gp_y * TILE_SIZE + ly; + int y_bottom = rows - 1 - y_top; + + #pragma unroll + for (int i = 0; i < TILE_SIZE; i += BLOCK_ROWS) + { + if (x < cols && y_top + i < half_rows) + { + int curr_y_top = y_top + i; + int curr_y_bottom = rows - 1 - curr_y_top; + + T val_top = loadpix(srcptr + mad24(curr_y_top, src_step, mad24(x, TSIZE, src_offset))); + T val_bottom = loadpix(srcptr + mad24(curr_y_bottom, src_step, mad24(x, TSIZE, src_offset))); + + tile_top[mad24(ly + i, LDS_STEP, lx)] = val_top; + tile_bottom[mad24(ly + i, LDS_STEP, lx)] = val_bottom; + } + } + barrier(CLK_LOCAL_MEM_FENCE); + + #pragma unroll + for (int i = 0; i < TILE_SIZE; i += BLOCK_ROWS) + { + if (x < cols && y_top + i < half_rows) + { + int curr_y_top = y_top + i; + int curr_y_bottom = rows - 1 - curr_y_top; + + storepix(tile_bottom[mad24(ly + i, LDS_STEP, lx)], + srcptr + mad24(curr_y_top, src_step, mad24(x, TSIZE, src_offset))); + + if (curr_y_top != curr_y_bottom) + { + storepix(tile_top[mad24(ly + i, LDS_STEP, lx)], + srcptr + mad24(curr_y_bottom, src_step, mad24(x, TSIZE, src_offset))); + } + } + } +} + +__kernel void arithm_flip_rows_cols_inplace(__global uchar * srcptr, int src_step, int src_offset, + int rows, int cols) +{ + int gp_x = get_group_id(0); + int gp_y = get_group_id(1); + int lx = get_local_id(0); + int ly = get_local_id(1); + + __local T tile_first[TILE_SIZE * LDS_STEP]; + __local T tile_second[TILE_SIZE * LDS_STEP]; + + int total_pixels = rows * cols; + int half_pixels = (total_pixels + 1) / 2; + + int x_first = gp_x * TILE_SIZE + lx; + int y_first = gp_y * TILE_SIZE + ly; + + int x_second = cols - 1 - x_first; + int y_second = rows - 1 - y_first; + + #pragma unroll + for (int i = 0; i < TILE_SIZE; i += BLOCK_ROWS) + { + int curr_y_first = y_first + i; + int curr_y_second = rows - 1 - curr_y_first; + int linear_idx = curr_y_first * cols + x_first; + + if (x_first < cols && curr_y_first < rows && linear_idx < half_pixels) + { + T val_first = loadpix(srcptr + mad24(curr_y_first, src_step, mad24(x_first, TSIZE, src_offset))); + T val_second = loadpix(srcptr + mad24(curr_y_second, src_step, mad24(x_second, TSIZE, src_offset))); + +#if kercn == 2 +#if cn == 1 + val_first = val_first.s10; + val_second = val_second.s10; +#endif +#elif kercn == 4 +#if cn == 1 + val_first = val_first.s3210; + val_second = val_second.s3210; +#elif cn == 2 + val_first = val_first.s2301; + val_second = val_second.s2301; +#endif +#endif + tile_first[mad24(ly + i, LDS_STEP, lx)] = val_first; + tile_second[mad24(ly + i, LDS_STEP, lx)] = val_second; + } + } + barrier(CLK_LOCAL_MEM_FENCE); + + #pragma unroll + for (int i = 0; i < TILE_SIZE; i += BLOCK_ROWS) + { + int curr_y_first = y_first + i; + int curr_y_second = rows - 1 - curr_y_first; + int linear_idx = curr_y_first * cols + x_first; + + if (x_first < cols && curr_y_first < rows && linear_idx < half_pixels) + { + storepix(tile_second[mad24(ly + i, LDS_STEP, lx)], + srcptr + mad24(curr_y_first, src_step, mad24(x_first, TSIZE, src_offset))); + + if (linear_idx != total_pixels - 1 - linear_idx) + { + storepix(tile_first[mad24(ly + i, LDS_STEP, lx)], + srcptr + mad24(curr_y_second, src_step, mad24(x_second, TSIZE, src_offset))); + } + } + } +} + +__kernel void arithm_flip_cols_inplace(__global uchar * srcptr, int src_step, int src_offset, + int rows, int cols) +{ + int gp_x = get_group_id(0); + int gp_y = get_group_id(1); + int lx = get_local_id(0); + int ly = get_local_id(1); + + __local T tile_left[TILE_SIZE * LDS_STEP]; + __local T tile_right[TILE_SIZE * LDS_STEP]; + + int half_cols = (cols + 1) / 2; + + int x_left = gp_x * TILE_SIZE + lx; + int y = gp_y * TILE_SIZE + ly; + int x_right = cols - 1 - x_left; + + #pragma unroll + for (int i = 0; i < TILE_SIZE; i += BLOCK_ROWS) + { + if (y + i < rows && x_left < half_cols) + { + T val_left = loadpix(srcptr + mad24(y + i, src_step, mad24(x_left, TSIZE, src_offset))); + T val_right = loadpix(srcptr + mad24(y + i, src_step, mad24(x_right, TSIZE, src_offset))); + +#if kercn == 2 +#if cn == 1 + val_left = val_left.s10; + val_right = val_right.s10; +#endif +#elif kercn == 4 +#if cn == 1 + val_left = val_left.s3210; + val_right = val_right.s3210; +#elif cn == 2 + val_left = val_left.s2301; + val_right = val_right.s2301; +#endif +#endif + tile_left[mad24(ly + i, LDS_STEP, lx)] = val_left; + tile_right[mad24(ly + i, LDS_STEP, lx)] = val_right; + } + } + barrier(CLK_LOCAL_MEM_FENCE); + + #pragma unroll + for (int i = 0; i < TILE_SIZE; i += BLOCK_ROWS) + { + if (y + i < rows && x_left < half_cols) + { + storepix(tile_right[mad24(ly + i, LDS_STEP, lx)], + srcptr + mad24(y + i, src_step, mad24(x_left, TSIZE, src_offset))); + + if (x_left != x_right) + { + storepix(tile_left[mad24(ly + i, LDS_STEP, lx)], + srcptr + mad24(y + i, src_step, mad24(x_right, TSIZE, src_offset))); + } + } + } +} + +#endif // INPLACE \ No newline at end of file From 9f101a126f1c7b0f71f56491ce751e009fce8fff Mon Sep 17 00:00:00 2001 From: 4ekmah Date: Mon, 20 Apr 2026 17:22:19 +0300 Subject: [PATCH 16/19] Merge pull request #28802 from 4ekmah:pyr_ecc Multiscale ECC #28802 OpenCV Extra: https://github.com/opencv/opencv_extra/pull/1338 ### 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 --- .../video/include/opencv2/video/tracking.hpp | 82 ++ modules/video/perf/perf_ecc.cpp | 60 ++ modules/video/src/eccms.cpp | 885 ++++++++++++++++++ modules/video/test/test_ecc.cpp | 368 ++++---- 4 files changed, 1211 insertions(+), 184 deletions(-) create mode 100644 modules/video/src/eccms.cpp diff --git a/modules/video/include/opencv2/video/tracking.hpp b/modules/video/include/opencv2/video/tracking.hpp index 2ed27f505d..f62660494c 100644 --- a/modules/video/include/opencv2/video/tracking.hpp +++ b/modules/video/include/opencv2/video/tracking.hpp @@ -419,6 +419,88 @@ CV_EXPORTS_W double findTransformECCWithMask( InputArray templateImage, TermCriteria criteria = TermCriteria(TermCriteria::COUNT + TermCriteria::EPS, 50, 1e-6), int gaussFiltSize = 5 ); +/** @brief struct ECCParameters is used by findTransformECCMultiScale + +@param motionType parameter, specifying the type of motion: + - **MOTION_TRANSLATION** sets a translational motion model; warpMatrix is \f$2\times 3\f$ with + the first \f$2\times 2\f$ part being the unity matrix and the rest two parameters being + estimated. + - **MOTION_EUCLIDEAN** sets a Euclidean (rigid) transformation as motion model; three + parameters are estimated; warpMatrix is \f$2\times 3\f$. + - **MOTION_AFFINE** sets an affine motion model (DEFAULT); six parameters are estimated; + warpMatrix is \f$2\times 3\f$. + - **MOTION_HOMOGRAPHY** sets a homography as a motion model; eight parameters are + estimated;\`warpMatrix\` is \f$3\times 3\f$. +@param criteria parameter, specifying the termination criteria of the ECC algorithm; +criteria.epsilon defines the threshold of the increment in the correlation coefficient between two +iterations (a negative criteria.epsilon makes criteria.maxcount the only termination criterion). +Default values are shown in the declaration above. +@param itersPerLevel Criterion extension: distribution of iterations limit over pyramid levels. +Can be empty, in this case, this algorithm will use criteria.maxCount on each level. +@param gaussFiltSize An optional value indicating size of gaussian blur filter; (DEFAULT: 5) +@param nlevels An optional value indicating amount of levels in the pyramid; (DEFAULT: 4) +@param interpolation Type of warp interpolation. Possible values are INTER_NEAREST and INTER_LINEAR. +Affects accuracy, especially when motionType == MOTION_TRANSLATION. (DEFAULT: INTER_LINEAR) + */ +struct CV_EXPORTS_W_SIMPLE ECCParameters +{ + CV_WRAP ECCParameters() {} + CV_PROP_RW int motionType = MOTION_AFFINE; + CV_PROP_RW cv::TermCriteria criteria = TermCriteria(TermCriteria::COUNT + TermCriteria::EPS, 50, 1e-6); + CV_PROP_RW std::vector itersPerLevel = std::vector(); + CV_PROP_RW int gaussFiltSize = 5; + CV_PROP_RW int nlevels = 4; + CV_PROP_RW int interpolation = INTER_LINEAR; +}; + +/** @brief Finds the geometric transform (warp) between two images in terms of the ECC criterion @cite EP08. Uses pyramids. + +@param reference Single channel reference image; CV_8U, CV_16U, CV_32F, CV_64F type. +@param sample sample image which should be warped with the final warpMatrix in +order to provide an image similar to reference, same type as reference. +@param warpMatrix floating-point \f$2\times 3\f$ or \f$3\times 3\f$ mapping matrix (warp). +@param eccParams List of the algorithm parameters. See ECCParameters for details. +@param referenceMask An optional single channel mask to indicate valid values of reference. +@param sampleMask An optional single channel mask to indicate valid values of sample. + +The function estimates the optimum transformation (warpMatrix) with respect to ECC criterion +(@cite EP08), that is + +\f[\texttt{warpMatrix} = \arg\max_{W} \texttt{ECC}(\texttt{templateImage}(x,y),\texttt{inputImage}(x',y'))\f] + +where + +\f[\begin{bmatrix} x' \\ y' \end{bmatrix} = W \cdot \begin{bmatrix} x \\ y \\ 1 \end{bmatrix}\f] + +(the equation holds with homogeneous coordinates for homography). It returns the final enhanced +correlation coefficient, that is the correlation coefficient between the template image and the +final warped input image. When a \f$3\times 3\f$ matrix is given with motionType =0, 1 or 2, the third +row is ignored. + +Unlike findHomography and estimateRigidTransform, the function findTransformECCMultiScale implements +an area-based alignment that builds on intensity similarities. In essence, the function updates the +initial transformation that roughly aligns the images. If this information is missing, the identity +warp (unity matrix) is used as an initialization. Note that if images undergo strong +displacements/rotations, an initial transformation that roughly aligns the images is necessary +(e.g., a simple euclidean/similarity transform that allows for the images showing the same image +content approximately). Use inverse warping in the second image to take an image close to the first +one, i.e. use the flag WARP_INVERSE_MAP with warpAffine or warpPerspective. See also the OpenCV +sample image_alignment.cpp that demonstrates the use of the function. Note that the function throws +an exception if algorithm does not converges. +Unlike findTransformECC, the findTransformECCMultiScale uses pyramids, making function more stable +and able to handle correctly more sophisticated cases. + +@sa +computeECC, estimateAffine2D, estimateAffinePartial2D, findHomography +*/ + +CV_EXPORTS_W double findTransformECCMultiScale(InputArray reference, + InputArray sample, + InputOutputArray warpMatrix, + const ECCParameters& eccParams = ECCParameters(), + InputArray referenceMask = noArray(), + InputArray sampleMask = noArray()); + /** @example samples/cpp/kalman.cpp An example using the standard Kalman filter */ diff --git a/modules/video/perf/perf_ecc.cpp b/modules/video/perf/perf_ecc.cpp index f81f64d8a6..87f1fbf32f 100644 --- a/modules/video/perf/perf_ecc.cpp +++ b/modules/video/perf/perf_ecc.cpp @@ -5,9 +5,14 @@ using namespace perf; CV_ENUM(MotionType, MOTION_TRANSLATION, MOTION_EUCLIDEAN, MOTION_AFFINE, MOTION_HOMOGRAPHY) CV_ENUM(ReadFlag, IMREAD_GRAYSCALE, IMREAD_COLOR) +CV_ENUM(MultiScaleFlag, false, true) typedef std::tuple TestParams; +typedef std::tuple TestParamsMS; typedef perf::TestBaseWithParam ECCPerfTest; +typedef perf::TestBaseWithParam ECCPerfTestMS; + +typedef std::tuple TestParams; PERF_TEST_P(ECCPerfTest, findTransformECC, testing::Combine(testing::Values(MOTION_TRANSLATION, MOTION_EUCLIDEAN, MOTION_AFFINE, MOTION_HOMOGRAPHY), @@ -67,4 +72,59 @@ PERF_TEST_P(ECCPerfTest, findTransformECC, } } +PERF_TEST_P(ECCPerfTestMS, findTransformECCMultiScale, + testing::Combine(testing::Values(MOTION_TRANSLATION, MOTION_EUCLIDEAN, MOTION_AFFINE, MOTION_HOMOGRAPHY), + testing::Values(false, true))) { + int transform_type = get<0>(GetParam()); + bool multiscaleFlag = get<1>(GetParam()); + + Mat img = imread(getDataPath("cv/shared/3MP.png"), IMREAD_GRAYSCALE); + Mat templateImage; + + Mat warpMat; + Mat warpGround; + + double angle; + switch (transform_type) { + case MOTION_TRANSLATION: + warpGround = (Mat_(2, 3) << 1.f, 0.f, 7.234f, 0.f, 1.f, 11.839f); + warpAffine(img, templateImage, warpGround, img.size(), INTER_LINEAR + WARP_INVERSE_MAP); + break; + case MOTION_EUCLIDEAN: + angle = CV_PI / 30; + + warpGround = (Mat_(2, 3) << (float)cos(angle), (float)-sin(angle), 12.123f, (float)sin(angle), + (float)cos(angle), 14.789f); + warpAffine(img, templateImage, warpGround, img.size(), INTER_LINEAR + WARP_INVERSE_MAP); + break; + case MOTION_AFFINE: + warpGround = (Mat_(2, 3) << 0.98f, 0.03f, 15.523f, -0.02f, 0.95f, 10.456f); + warpAffine(img, templateImage, warpGround, img.size(), INTER_LINEAR + WARP_INVERSE_MAP); + break; + case MOTION_HOMOGRAPHY: + warpGround = (Mat_(3, 3) << 0.98f, 0.03f, 15.523f, -0.02f, 0.95f, 10.456f, 0.0002f, 0.0003f, 1.f); + warpPerspective(img, templateImage, warpGround, img.size(), INTER_LINEAR + WARP_INVERSE_MAP); + break; + } + + TEST_CYCLE() { + if (transform_type < 3) + warpMat = Mat::eye(2, 3, CV_32F); + else + warpMat = Mat::eye(3, 3, CV_32F); + + if(multiscaleFlag) { + ECCParameters params; + params.criteria = cv::TermCriteria(TermCriteria::COUNT + TermCriteria::EPS, 5, -1); + params.motionType = transform_type; + params.itersPerLevel = {1, 2, 2, 2}; + findTransformECCMultiScale(templateImage, img, warpMat, params); + } + else { + findTransformECC(templateImage, img, warpMat, transform_type, + TermCriteria(TermCriteria::COUNT + TermCriteria::EPS, 5, -1)); + } + } + SANITY_CHECK_NOTHING(); +} } // namespace opencv_test diff --git a/modules/video/src/eccms.cpp b/modules/video/src/eccms.cpp new file mode 100644 index 0000000000..963763d5d3 --- /dev/null +++ b/modules/video/src/eccms.cpp @@ -0,0 +1,885 @@ +// 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 "precomp.hpp" + +/****************************************************************************************\ +* Image Alignment (ECC algorithm, pyramidal version) * +\****************************************************************************************/ + +namespace cv { +typedef std::vector MatPyramid; + +template struct MotionTraits {}; + +template<> struct MotionTraits { + enum { paramAmount = 2 }; + static inline void tailHandlerGetCoord(float& sx, + float& sy, + float& denominator, + int col, + float numeratorX0, + float numeratorY0, + float /*denominator0*/, + float /*a00*/, + float /*a10*/, + float /*a20*/) + { + denominator = 0; + sx = (numeratorX0 + col); + sy = numeratorY0; + } + template + static constexpr std::array fillJacobian(int /*col*/, int /*row*/, float/*sx*/, float/*sy*/, float fVal, + elemtype gx, elemtype gy, float /*a00*/, float /*a10*/, + float/*denominator*/) { +#define GX (fVal * gx) +#define GY (fVal * gy) + return std::array{GX, GY}; +#undef GX +#undef GY + } +}; + +template<> struct MotionTraits { + enum { paramAmount = 3 }; + static inline void tailHandlerGetCoord(float& sx, + float& sy, + float& denominator, + int col, + float numeratorX0, + float numeratorY0, + float /*denominator0*/, + float a00, + float a10, + float /*a20*/) + { + denominator = 0; + sx = (numeratorX0 + a00 * col); + sy = (numeratorY0 + a10 * col); + } + + template + static constexpr std::array fillJacobian(int col, int row, float/*sx*/, float/*sy*/, float fVal, + elemtype gx, elemtype gy, float a00, float a10, + float/*denominator*/) { +#define GX (fVal * gx) +#define GY (fVal * gy) +#define HATX (-col * a10 - row * a00) +#define HATY (col * a00 - row * a10) +#define GZ (GX * HATX + GY * HATY) + return std::array{GZ, GX, GY}; +#undef GX +#undef GY +#undef HATX +#undef HATY +#undef GZ + } +}; + +template<> struct MotionTraits { + enum { paramAmount = 6}; + static inline void tailHandlerGetCoord(float& sx, + float& sy, + float& denominator, + int col, + float numeratorX0, + float numeratorY0, + float /*denominator0*/, + float a00, + float a10, + float /*a20*/) + { + denominator = 0; + sx = (numeratorX0 + a00 * col); + sy = (numeratorY0 + a10 * col); + } + + template + static constexpr std::array fillJacobian(int col, int row, float/*sx*/, float/*sy*/, float fVal, + elemtype gx, elemtype gy, float /*a00*/, float /*a10*/, + float/*denominator*/) { +#define GX (fVal * gx) +#define GY (fVal * gy) + return std::array{GX * col, GY * col, GX * row, GY * row, GX, GY}; +#undef GX +#undef GY + } +}; + +template<> struct MotionTraits { + enum { paramAmount = 8}; + static inline void tailHandlerGetCoord(float& sx, + float& sy, + float& denominator, + int col, + float numeratorX0, + float numeratorY0, + float denominator0, + float a00, + float a10, + float a20) + { + denominator = 1.f / (col * a20 + denominator0); + sx = (numeratorX0 + a00 * col) * denominator; + sy = (numeratorY0 + a10 * col) * denominator; + } + + template + static constexpr std::array fillJacobian(int col, int row, float sx, float sy, float fVal, + elemtype gx, elemtype gy, float/*a00*/, float/*a10*/, + float denominator) { +#define GX (fVal * float(gx) * denominator) +#define GY (fVal * float(gy) * denominator) +#define GZ (-(GX * sx + GY * sy)) + return std::array{GX * col, GY * col, GZ * col, GX * row, GY * row, GZ * row, GX, GY}; +#undef GX +#undef GY +#undef GZ + } +}; + +inline void reinterpret(Mat& mat, int newdepth) { + mat.flags = (mat.flags & ~CV_MAT_DEPTH_MASK) | newdepth; +} + +template +class constexprForClass +{ +public: + static inline void execute(F&& fVal) { + constexprForClass::execute(std::forward(fVal)); + fVal(N-1); + } +}; + +template +class constexprForClass<0, F> +{ +public: + static inline void execute(F&&) {} +}; + +template +void constexprFor(F&& fVal) { + constexprForClass::execute(std::forward(fVal)); +} +template +class constexprForUpperTriangleClassOneRow +{ +public: + static inline void execute(F&& fVal) { + constexprForUpperTriangleClassOneRow::execute(std::forward(fVal)); + fVal(R, R + C - 1); + } +}; + +template +class constexprForUpperTriangleClassOneRow +{ +public: + static inline void execute(F&&) {} +}; + +template +class constexprForUpperTriangleClass +{ +public: + static inline void execute(F&& fVal) { + constexprForUpperTriangleClass::execute(std::forward(fVal)); + constexprForUpperTriangleClassOneRow::execute(std::forward(fVal)); + } +}; + +template +class constexprForUpperTriangleClass<0, D, F> +{ +public: + static inline void execute(F&&) {} +}; + +template +void constexprForUpperTriangle(F&& fVal) { + constexprForUpperTriangleClass::execute(std::forward(fVal)); +} + +template +constexpr int hessianRowStart(int row) { + return row == 0 ? 0 : (MotionTraits::paramAmount - row + 1 + hessianRowStart(row - 1)); +} + +template +static double imageHessianProjECC(const Mat& map, + const Mat& sampleWithGrad, + const Mat& ref, + double& sampSum, + double& sampSqSum, + double& refSum, + double& refSqSum, + int& nz, + Mat& hessian, + Mat& sampleProj, + Mat& refProj, + int deltaY, + int interpolation) { + static_assert(std::is_same::value, "imageHessianProjECC: f16 is not supported yet"); +#define HESSIAN_PARAMS (MotionTraits::paramAmount) + + CV_Assert(map.type() == CV_64F); + CV_Assert(interpolation == INTER_NEAREST || interpolation == INTER_LINEAR); + CV_Assert(hessian.type() == CV_64F && sampleProj.type() == CV_64F && refProj.type() == CV_64F); + if (sampleProj.size() != Size(1, HESSIAN_PARAMS) || refProj.size() != Size(1, HESSIAN_PARAMS)) { + CV_Error(Error::BadImageSize, format("imageHessianProjECC: Wrong sample projection/reference projection size. 1x%d expected", HESSIAN_PARAMS)); + } + if (hessian.size() != Size(HESSIAN_PARAMS, HESSIAN_PARAMS)) { + CV_Error(Error::BadImageSize, format("imageHessianProjECC: Wrong hessian size. %dx%d expected", HESSIAN_PARAMS, HESSIAN_PARAMS)); + } + if (!map.isContinuous()) { + CV_Error(Error::BadStep, "imageHessianProjECC: Map should be continuous"); + } + if (std::is_same::value) { + CV_Assert(sampleWithGrad.type() == CV_32FC4 && ref.type() == CV_32FC2); + } + + int hr = ref.rows; + int wr = ref.cols; + int hs = sampleWithGrad.rows; + int ws = sampleWithGrad.cols; + unsigned int ycond = hs - (INTER_LINEAR ? 1 : 0); + unsigned int xcond = ws - (INTER_LINEAR ? 1 : 0); + + hessian = Mat::zeros(hessian.size(), hessian.type()); + sampleProj = Mat::zeros(sampleProj.size(), sampleProj.type()); + refProj = Mat::zeros(refProj.size(), refProj.type()); + + const int MAX_STRIPES = 128; + int stripesAmount = std::min(MAX_STRIPES, hr / deltaY); + std::vector::paramAmount, MotionTraits::paramAmount> > hessPs(stripesAmount); + std::vector::paramAmount> > iprojs(stripesAmount); + std::vector::paramAmount> > tprojs(stripesAmount); + std::vector::paramAmount> > projSubs(stripesAmount); + std::vector correlations(stripesAmount, 0.); + + std::vector sampSums(stripesAmount, 0); + std::vector sampSqSums(stripesAmount, 0); + std::vector refSums(stripesAmount, 0); + std::vector refSqSums(stripesAmount, 0); + std::vector nzs(stripesAmount, 0); + std::vector sampMaskedSums(stripesAmount, 0); + std::vector refMaskedSums(stripesAmount, 0); + + double a00 = map.at(0, 0); + double a01 = map.at(0, 1); + double a02 = map.at(0, 2); + double a10 = map.at(1, 0); + double a11 = map.at(1, 1); + double a12 = map.at(1, 2); + double a20 = 0; + double a21 = 0; + double a22 = 0; + if (motionType == MOTION_HOMOGRAPHY) { + a20 = map.at(2, 0); + a21 = map.at(2, 1); + a22 = map.at(2, 2); + } + + const elemtype* samplePtr0 = sampleWithGrad.ptr(0); + + parallel_for_(Range(0, stripesAmount), [&](const Range& range) { + int stripeIdx = range.start; + int ystart = (hr * stripeIdx) / stripesAmount; + ystart = roundUp(ystart, deltaY); + int yend = (hr * (range.end)) / stripesAmount; + // we don't store intermediate jacobian; instead, we iteratively update Hessian, sampleProj and refProj + for (int y = ystart; y < yend; y += deltaY) { + const elemtype* refPtr = ref.ptr(y); + + std::array hessPcache{}; + std::array iprojCache{}; + std::array tprojCache{}; + std::array projSubCache{}; + + const float numeratorX0 = y * (float)a01 + (float)a02; + const float numeratorY0 = y * (float)a11 + (float)a12; + const float denominator0 = y * (float)a21 + (float)a22; + int x = 0; + for (; x < wr; x++) { //Tail handler + float sx, sy, denominator; + MotionTraits::tailHandlerGetCoord(sx, sy, denominator, x, numeratorX0, numeratorY0, + denominator0, (float)a00, (float)a10, (float)a20); + const unsigned int x0 = (interpolation == INTER_LINEAR) ? static_cast(std::floor(sx)) : saturate_cast(sx); + const unsigned int y0 = (interpolation == INTER_LINEAR) ? static_cast(std::floor(sy)) : saturate_cast(sy); + if(interpolation == INTER_LINEAR && (static_cast(x0 < xcond) & static_cast(y0 < ycond)) == 0) + continue; + if (interpolation == INTER_NEAREST && (static_cast(x0 < xcond) & static_cast(y0 < ycond)) == 0) + continue; + float sampleVal = 0; + float gx = 0; + float gy = 0; + float fVal = 0; + if(interpolation == INTER_LINEAR) { + const int x1 = x0 + 1; + const int y1 = y0 + 1; + + const float dx = sx - x0; + const float dy = sy - y0; + + const float p00_val = samplePtr0[4 * y0 * ws + 4 * x0]; + const float p01_val = samplePtr0[4 * y0 * ws + 4 * x1]; + const float p10_val = samplePtr0[4 * y1 * ws + 4 * x0]; + const float p11_val = samplePtr0[4 * y1 * ws + 4 * x1]; + const float p0_val = p00_val * (1.0f - dx) + p01_val * dx; + const float p1_val = p10_val * (1.0f - dx) + p11_val * dx; + + const float p00_gx = samplePtr0[4 * y0 * ws + 4 * x0 + 1]; + const float p01_gx = samplePtr0[4 * y0 * ws + 4 * x1 + 1]; + const float p10_gx = samplePtr0[4 * y1 * ws + 4 * x0 + 1]; + const float p11_gx = samplePtr0[4 * y1 * ws + 4 * x1 + 1]; + const float p0_gx = p00_gx * (1.0f - dx) + p01_gx * dx; + const float p1_gx = p10_gx * (1.0f - dx) + p11_gx * dx; + + const float p00_gy = samplePtr0[4 * y0 * ws + 4 * x0 + 2]; + const float p01_gy = samplePtr0[4 * y0 * ws + 4 * x1 + 2]; + const float p10_gy = samplePtr0[4 * y1 * ws + 4 * x0 + 2]; + const float p11_gy = samplePtr0[4 * y1 * ws + 4 * x1 + 2]; + const float p0_gy = p00_gy * (1.0f - dx) + p01_gy * dx; + const float p1_gy = p10_gy * (1.0f - dx) + p11_gy * dx; + + const float p00_mask = samplePtr0[4 * y0 * ws + 4 * x0 + 2] == 0.f ? 0.f : 1.f; + const float p01_mask = samplePtr0[4 * y0 * ws + 4 * x1 + 2] == 0.f ? 0.f : 1.f; + const float p10_mask = samplePtr0[4 * y1 * ws + 4 * x0 + 2] == 0.f ? 0.f : 1.f; + const float p11_mask = samplePtr0[4 * y1 * ws + 4 * x1 + 2] == 0.f ? 0.f : 1.f; + + sampleVal = p0_val * (1.0f - dy) + p1_val * dy; + gx = p0_gx * (1.0f - dy) + p1_gx * dy; + gy = p0_gy * (1.0f - dy) + p1_gy * dy; + fVal = p00_mask * p01_mask * p10_mask * p11_mask; + } + else { // if(interpolation == INTER_NEAREST) + const elemtype* samplePtr = samplePtr0 + y0 * (ws * 4) + x0 * 4; + sampleVal = samplePtr[0]; + gx = samplePtr[1]; + gy = samplePtr[2]; + fVal = float(samplePtr[3]) == 0.f ? 0.f : 1.f; + } + + float refVal = refPtr[2 * x]; + fVal *= float(refPtr[2 * x + 1]) == 0.f ? 0.f : 1.f; + sampleVal *= fVal; + refVal *= fVal; + sampSums[stripeIdx] += sampleVal; + sampSqSums[stripeIdx] += sampleVal * sampleVal; + refSums[stripeIdx] += refVal; + refSqSums[stripeIdx] += refVal * refVal; + nzs[stripeIdx] += (int)fVal; + sampMaskedSums[stripeIdx] += sampleVal; + refMaskedSums[stripeIdx] += refVal; + std::array jac = MotionTraits::fillJacobian(x, y, sx, sy, + fVal, gx, + gy, (float)a00, + (float)a10, denominator); + constexprForUpperTriangle([&](int row_i, int col_i) { + hessPcache[hessianRowStart(row_i) + (col_i - row_i)] += jac[row_i] * jac[col_i]; + }); + constexprFor([&](int elem) { + iprojCache[elem] += jac[elem] * sampleVal; + tprojCache[elem] += jac[elem] * refVal; + projSubCache[elem] += jac[elem] * fVal; + }); + correlations[stripeIdx] += sampleVal * refVal; + } + + constexprForUpperTriangle([&](int row, int col) { + hessPs[stripeIdx](row, col) += hessPcache[hessianRowStart(row) + (col - row)]; + }); + constexprFor([&](int elem) { + iprojs[stripeIdx][elem] += iprojCache[elem]; + tprojs[stripeIdx][elem] += tprojCache[elem]; + projSubs[stripeIdx][elem] += projSubCache[elem]; + }); + } + }); + double sampMaskedSum = 0; + double refMaskedSum = 0; + double correlation = 0; + sampSum = sampSqSum = refSum = refSqSum = nz = 0; + + for (int stripeIdx = 0; stripeIdx < stripesAmount; stripeIdx++) { + correlation += correlations[stripeIdx]; + sampSum += sampSums[stripeIdx]; + sampSqSum += sampSqSums[stripeIdx]; + refSum += refSums[stripeIdx]; + refSqSum += refSqSums[stripeIdx]; + sampMaskedSum += sampMaskedSums[stripeIdx]; + refMaskedSum += refMaskedSums[stripeIdx]; + nz += nzs[stripeIdx]; + } + double scale = nz == 0 ? 0. : 1. / nz; + double sampMean = sampSum * scale; + double refMean = refSum * scale; + correlation += nz * sampMean * refMean - sampMaskedSum * refMean - refMaskedSum * sampMean; + double* hessPtr = hessian.ptr(0); + double* sampleProjPtr = sampleProj.ptr(0); + double* refProjPtr = refProj.ptr(0); + for (int stripeIdx = 0; stripeIdx < stripesAmount; stripeIdx++) { + for (int hessNum = 0; hessNum < HESSIAN_PARAMS * HESSIAN_PARAMS; hessNum++) { + hessPtr[hessNum] += hessPs[stripeIdx].val[hessNum]; + } + for (int projNum = 0; projNum < HESSIAN_PARAMS; projNum++) { + sampleProjPtr[projNum] += iprojs[stripeIdx][projNum] - projSubs[stripeIdx][projNum] * sampMean; + refProjPtr[projNum] += tprojs[stripeIdx][projNum] - projSubs[stripeIdx][projNum] * refMean; + } + } + constexprForUpperTriangle([&](int row, int col) { + hessPtr[col * HESSIAN_PARAMS + row] = hessPtr[row * HESSIAN_PARAMS + col]; + }); + return correlation; +#undef HESSIAN_PARAMS +} + +static void updateWarpingMatrixECC(Mat& map_matrix, const Mat& update, const int motionType) { + CV_Assert(map_matrix.type() == CV_64FC1); + CV_Assert(update.type() == CV_64FC1); + + CV_Assert(motionType == MOTION_TRANSLATION || motionType == MOTION_EUCLIDEAN || motionType == MOTION_AFFINE || + motionType == MOTION_HOMOGRAPHY); + + if (motionType == MOTION_HOMOGRAPHY) + CV_Assert(map_matrix.rows == 3 && update.rows == 8); + else if (motionType == MOTION_AFFINE) + CV_Assert(map_matrix.rows == 2 && update.rows == 6); + else if (motionType == MOTION_EUCLIDEAN) + CV_Assert(map_matrix.rows == 2 && update.rows == 3); + else + CV_Assert(map_matrix.rows == 2 && update.rows == 2); + + CV_Assert(update.cols == 1); + + CV_Assert(map_matrix.isContinuous()); + CV_Assert(update.isContinuous()); + + double* mapPtr = map_matrix.ptr(0); + const double* updatePtr = update.ptr(0); + + if (motionType == MOTION_TRANSLATION) { + mapPtr[2] += updatePtr[0]; + mapPtr[5] += updatePtr[1]; + } + if (motionType == MOTION_AFFINE) { + mapPtr[0] += updatePtr[0]; + mapPtr[3] += updatePtr[1]; + mapPtr[1] += updatePtr[2]; + mapPtr[4] += updatePtr[3]; + mapPtr[2] += updatePtr[4]; + mapPtr[5] += updatePtr[5]; + } + if (motionType == MOTION_HOMOGRAPHY) { + mapPtr[0] += updatePtr[0]; + mapPtr[3] += updatePtr[1]; + mapPtr[6] += updatePtr[2]; + mapPtr[1] += updatePtr[3]; + mapPtr[4] += updatePtr[4]; + mapPtr[7] += updatePtr[5]; + mapPtr[2] += updatePtr[6]; + mapPtr[5] += updatePtr[7]; + } + if (motionType == MOTION_EUCLIDEAN) { + double new_theta = updatePtr[0]; + new_theta += asin(mapPtr[3]); + + mapPtr[2] += updatePtr[1]; + mapPtr[5] += updatePtr[2]; + mapPtr[0] = mapPtr[4] = cos(new_theta); + mapPtr[3] = sin(new_theta); + mapPtr[1] = -mapPtr[3]; + } +} + +static void optimizeECC(Mat& sampleWithGrad, + const Mat& reference, + Mat& map, + int motionType, + double* rho, + double* lastRho, + int deltaY, + int nparams, + int interpolation) { + CV_Assert(interpolation == INTER_NEAREST || interpolation == INTER_LINEAR); + + // warp-back portion of the inputImage and gradients to the coordinate space of the referenceFloat + double correlation = 0; + + // matrices needed for solving linear equation system for maximizing ECC + Mat hessian = Mat(nparams, nparams, CV_64F); + Mat hessianInv = Mat(nparams, nparams, CV_64F); + Mat sampleProjection = Mat(nparams, 1, CV_64F); + Mat referenceProjection = Mat(nparams, 1, CV_64F); + Mat sampleProjectionHessian = Mat(nparams, 1, CV_64F); + Mat errorProjection = Mat(nparams, 1, CV_64F); + Mat deltaP = Mat(nparams, 1, CV_64F); + + double sampSum; + double sampSqSum; + double referenceSum; + double referenceSqSum; + int nz; + + { // if(imageWithGrad.type() == CV_32FC4) + if (motionType == MOTION_TRANSLATION) { + correlation = imageHessianProjECC(map, + sampleWithGrad, + reference, + sampSum, + sampSqSum, + referenceSum, + referenceSqSum, + nz, + hessian, + sampleProjection, + referenceProjection, + deltaY, + interpolation); + } else if (motionType == MOTION_EUCLIDEAN) { + correlation = imageHessianProjECC(map, + sampleWithGrad, + reference, + sampSum, + sampSqSum, + referenceSum, + referenceSqSum, + nz, + hessian, + sampleProjection, + referenceProjection, + deltaY, + interpolation); + } else if (motionType == MOTION_AFFINE) { + correlation = imageHessianProjECC(map, + sampleWithGrad, + reference, + sampSum, + sampSqSum, + referenceSum, + referenceSqSum, + nz, + hessian, + sampleProjection, + referenceProjection, + deltaY, + interpolation); + } else { + correlation = imageHessianProjECC(map, + sampleWithGrad, + reference, + sampSum, + sampSqSum, + referenceSum, + referenceSqSum, + nz, + hessian, + sampleProjection, + referenceProjection, + deltaY, + interpolation); + } + } + double scale = nz == 0 ? 0. : 1. / nz; + double sampMean = sampSum * scale; + double refMean = referenceSum * scale; + double sampStd = std::sqrt(std::max(sampSqSum * scale - sampMean * sampMean, 0.)); + double refStd = std::sqrt(std::max(referenceSqSum * scale - refMean * refMean, 0.)); + + // inverse of Hessian + hessianInv = hessian.inv(); + // calculate enhanced correlation coefficient (ECC)->rho + *lastRho = *rho; + double refNorm = std::sqrt(nz * refStd * refStd); + double sampNorm = std::sqrt(nz * sampStd * sampStd); + + *rho = correlation / (sampNorm * refNorm); + if ((bool)cvIsNaN(*rho)) { + CV_Error(Error::StsNoConv, "NaN encountered."); + } + + // calculate the parameter lambda to account for illumination variation + sampleProjectionHessian = hessianInv * sampleProjection; + const double lambdaN = (sampNorm * sampNorm) - sampleProjection.dot(sampleProjectionHessian); + const double lambdaD = correlation - referenceProjection.dot(sampleProjectionHessian); + + if (lambdaD <= 0.0) { + CV_Error(Error::StsNoConv, "The algorithm stopped before its convergence. The correlation is going to be minimized. " + "Images may be uncorrelated or non-overlapped"); + } + const double lambda = (lambdaN / lambdaD); + + // estimate the update step delta_p + errorProjection = lambda * referenceProjection - sampleProjection; + gemm(hessianInv, errorProjection, 1., noArray(), 0., deltaP); + + // update warping matrix + updateWarpingMatrixECC(map, deltaP, motionType); +} + +static Mat prepareGradients(const Mat& sample) { + CV_Assert(sample.type() == CV_32FC2 || sample.type() == CV_16FC2); + + const int ws = sample.cols; + const int hs = sample.rows; + + Mat sampleWithGrad; + int ntasks = std::min(4, hs); + + { + sampleWithGrad = Mat(hs, ws, CV_32FC4); + float* dstPtr = sampleWithGrad.ptr(); + parallel_for_(Range(0, ntasks), [&](const Range& range) { + int rowstart = range.start * hs / ntasks; + int rowend = range.end * hs / ntasks; + for (int row = rowstart; row < rowend; row++) { + const float* sampleCurLine = sample.ptr(row); + const float* samplePrevLine = sample.ptr(std::max(row - 1, 0)); + const float* sampleNextLine = sample.ptr(std::min(row + 1, hs - 1)); + float gradDivY = (row > 0 && row + 1 < hs) ? 0.5f : 0.25f; + int col = 0; + for (; col < ws; col++) { + int prevCol = std::max(col - 1, 0); + int nextCol = std::min(col + 1, ws - 1); + float gradDivX = (col > 0 && col + 1 < ws) ? 0.5f : 0.25f; + dstPtr[row * ws * 4 + col * 4] = sampleCurLine[2 * col]; + dstPtr[row * ws * 4 + col * 4 + 1] = + gradDivX * (sampleCurLine[2 * nextCol] - sampleCurLine[2 * prevCol]); + dstPtr[row * ws * 4 + col * 4 + 2] = gradDivY * (sampleNextLine[2 * col] - samplePrevLine[2 * col]); + dstPtr[row * ws * 4 + col * 4 + 3] = sampleCurLine[2 * col + 1]; + } + } + }, ntasks); + } + return sampleWithGrad; +} + +static void buildPyramidECC(InputArray inputImage, + MatPyramid& imgPyramid, + InputArray& mask, + MatPyramid& maskPyramid, + int nlevels) { + imgPyramid.resize(nlevels); + inputImage.getMat().convertTo(imgPyramid[0], CV_8UC1); + maskPyramid.resize(nlevels); + if (!mask.empty()) { + mask.getMat().convertTo(maskPyramid[0], CV_8UC1); + } + for (int pyrLevel = 0; pyrLevel < nlevels - 1; ++pyrLevel) { + Size size = Size((imgPyramid[pyrLevel].cols + 1) / 2, (imgPyramid[pyrLevel].rows + 1) / 2); + pyrDown(imgPyramid[pyrLevel], imgPyramid[pyrLevel + 1], size); + if (!mask.empty()) { + pyrDown(maskPyramid[pyrLevel], maskPyramid[pyrLevel + 1], size); + threshold(maskPyramid[pyrLevel + 1], maskPyramid[pyrLevel + 1], 254, 0xff, THRESH_BINARY); + } + } +} + +static Mat spliceWithMask(const Mat& image, const Mat& mask) { + CV_Assert(image.type() == CV_32F && (mask.empty() || mask.type() == CV_8U)); + if (!mask.empty() && image.size() != mask.size()) { + CV_Error(Error::BadImageSize, "spliceWithMask: Mask and image have to be of same size."); + } + const int hs = image.rows; + const int ws = image.cols; + + Mat result; + int ntasks = std::min(4, hs); + { + union conv_ { + uint32_t valU; + float val; + conv_() : valU(0xffffffff) {} + } conv; + result = Mat(hs, ws, CV_32FC2); + parallel_for_(Range(0, ntasks), [&](const Range& range) { + int rowstart = range.start * hs / ntasks; + int rowend = range.end * hs / ntasks; + for (int row = rowstart; row < rowend; row++) { + float* dstPtr = result.ptr(row); + const float* srcPtr = image.ptr(row); + const uint8_t* maskPtr = !mask.empty() ? mask.ptr(row) : nullptr; + int col = 0; + for (; col < ws; col++) { + dstPtr[col * 2] = srcPtr[col]; + dstPtr[col * 2 + 1] = (!maskPtr || maskPtr[col]) ? conv.val : 0; + } + } + }, ntasks); + } + return result; +} + +static void scaleWarpMatrix(Mat& warpMatrix, float scale) { + if (warpMatrix.rows == 3) { + Mat invertScaleMat = Mat(3, 3, CV_64F, 0.f); + invertScaleMat.at(0, 0) = 1.f / scale; + invertScaleMat.at(1, 1) = 1.f / scale; + invertScaleMat.at(2, 2) = 1.f; + Mat scaleMatrix = invertScaleMat.clone(); + scaleMatrix.at(0, 0) = scale; + scaleMatrix.at(1, 1) = scale; + gemm(warpMatrix, invertScaleMat, 1., noArray(), 0., warpMatrix); + gemm(scaleMatrix, warpMatrix, 1., noArray(), 0., warpMatrix); + // Normalization, internal algorithms assumes, that a22 = 1.0f + for (int mel = 0; mel < 8; mel++) { + (reinterpret_cast(warpMatrix.data))[mel] /= + (reinterpret_cast(warpMatrix.data))[8]; + } + (reinterpret_cast(warpMatrix.data))[8] = 1.f; + } else { + warpMatrix.at(0, 2) *= scale; + warpMatrix.at(1, 2) *= scale; + } +} + +static void checkParams(const MatPyramid& referencePyramid, + const MatPyramid& samplePyramid, + Mat& map, + std::vector& itersPerLevel, + const ECCParameters& eccParams) { + if (itersPerLevel.empty()) { + itersPerLevel.resize(eccParams.nlevels, eccParams.criteria.maxCount); + } + CV_Assert(eccParams.interpolation == INTER_NEAREST || eccParams.interpolation == INTER_LINEAR); + CV_Assert(static_cast(itersPerLevel.size()) == eccParams.nlevels); + for (const auto& lvl : referencePyramid) { + CV_Assert(!lvl.empty() && lvl.type() == referencePyramid[0].type()); + } + CV_Assert(!samplePyramid.empty()); + for (const auto& lvl : samplePyramid) { + CV_Assert(!lvl.empty() && lvl.type() == samplePyramid[0].type()); + } + CV_Assert(samplePyramid.size() == referencePyramid.size() && samplePyramid.size() == itersPerLevel.size()); + CV_Assert(referencePyramid.back().rows > 1 && referencePyramid.back().cols > 1 && + samplePyramid.back().rows > 1 && samplePyramid.back().cols > 1); + // If the user passed an un-initialized warpMatrix, initialize to identity + if (referencePyramid[0].type() != CV_32FC2 && referencePyramid[0].type() != CV_16FC2) { + CV_Error(Error::StsError, "Reference pyramid have to be prepared via prepareReferencePyramid function"); + } + // accept only 1-channel images + CV_Assert(samplePyramid[0].type() == CV_32FC2 || samplePyramid[0].type() != CV_16FC2); + CV_Assert(map.type() == CV_64FC1); + if (map.cols != 3 || (map.rows != 2 && map.rows != 3)) { + CV_Error(Error::BadImageSize, "warpMatrix has incorrect size"); + } + + if (eccParams.motionType != MOTION_TRANSLATION && eccParams.motionType != MOTION_EUCLIDEAN && + eccParams.motionType != MOTION_AFFINE && eccParams.motionType != MOTION_HOMOGRAPHY) { + CV_Error(Error::StsError, "Incorrect motion type"); + } + + if (eccParams.motionType == MOTION_HOMOGRAPHY && map.rows != 3) { + CV_Error(Error::BadImageSize, "warpMatrix has incorrect size"); + } + + if (!((bool)(eccParams.criteria.type & TermCriteria::COUNT) || (bool)(eccParams.criteria.type & TermCriteria::EPS))) { + CV_Error(Error::StsError, "Incorrect stop eccParams.criteria"); + } +} + +static MatPyramid prepareECCPyramid(InputArray image, + InputArray imageMask, // Can be empty + int gaussFiltSize, + int nlevels) { + MatPyramid imagePyramid, maskPyramid; + buildPyramidECC(image, imagePyramid, imageMask, maskPyramid, nlevels); + for (int lvl = 0; lvl < nlevels; lvl++) { + Mat imgFloat; + imagePyramid[lvl].convertTo(imgFloat, CV_32F, 1. / 255.); + if (gaussFiltSize != 0) { + GaussianBlur(imgFloat, imgFloat, Size(gaussFiltSize, gaussFiltSize), 0, 0); + } + imagePyramid[lvl] = spliceWithMask( + imgFloat, + (static_cast(maskPyramid.size()) > lvl && !maskPyramid[lvl].empty()) ? maskPyramid[lvl] : Mat()); + } + return imagePyramid; +} + +double findTransformECCMultiScale(InputArray reference, + InputArray sample, + InputOutputArray warpMatrixA, + const ECCParameters& eccParams, + InputArray referenceMask, + InputArray sampleMask) { + MatPyramid referencePyramid = prepareECCPyramid(reference, referenceMask, eccParams.gaussFiltSize, eccParams.nlevels); + MatPyramid samplePyramid = prepareECCPyramid(sample, sampleMask, eccParams.gaussFiltSize, eccParams.nlevels); + Mat& warpMatrix = warpMatrixA.getMatRef(); + std::vector itersPerLevelCopy = eccParams.itersPerLevel; + // If the user passed an un-initialized warpMatrix, initialize to identity + if (warpMatrix.empty()) + { + int rowCount = eccParams.motionType == MOTION_HOMOGRAPHY ? 3 : 2; + warpMatrix = Mat::eye(rowCount, 3, CV_64FC1); + } + int warpMatrixType = warpMatrix.type(); + if (warpMatrixType != CV_64FC1) + { + warpMatrix.convertTo(warpMatrix, CV_64FC1); + } + + checkParams(referencePyramid, + samplePyramid, + warpMatrix, + itersPerLevelCopy, + eccParams); + + int nparams = 0; + switch (eccParams.motionType) { + case MOTION_TRANSLATION: + nparams = MotionTraits::paramAmount; + break; + case MOTION_EUCLIDEAN: + nparams = MotionTraits::paramAmount; + break; + case MOTION_AFFINE: + nparams = MotionTraits::paramAmount; + break; + case MOTION_HOMOGRAPHY: + nparams = MotionTraits::paramAmount; + break; + default: + CV_Error(Error::StsBadArg, "Incorrect motion type"); + } + + const std::vector numberOfIterations = ((eccParams.criteria.type & TermCriteria::COUNT) != 0) + ? itersPerLevelCopy + : std::vector(eccParams.nlevels, 200); + const double terminationEPS = (bool)(eccParams.criteria.type & TermCriteria::EPS) ? eccParams.criteria.epsilon : -1; + + // Scale warp matrix multiple times to lower pyramid level + for (int pyrLevel = 0; pyrLevel < eccParams.nlevels - 1; pyrLevel++) { + scaleWarpMatrix(warpMatrix, 0.5); + } + double rho = -1; + for (int pyrLevel = eccParams.nlevels - 1; pyrLevel >= 0; --pyrLevel) { + const int hr = referencePyramid[pyrLevel].rows; + + Mat sampleWithGrad = prepareGradients(samplePyramid[pyrLevel]); + + const int LOW_SIZE = 200; + int deltaY = hr < LOW_SIZE ? 1 : 2; + + // iteratively update mapMatrix + double lastRho = -terminationEPS; + for (int i = 1; (i <= numberOfIterations[pyrLevel]) && (fabs(rho - lastRho) >= terminationEPS); i++) { + optimizeECC(sampleWithGrad, referencePyramid[pyrLevel], warpMatrix, eccParams.motionType, &rho, &lastRho, deltaY, nparams, eccParams.interpolation); + } + if (pyrLevel > 0) { + scaleWarpMatrix(warpMatrix, 2); + } + } + if(warpMatrixType != CV_64FC1) + { + warpMatrix.convertTo(warpMatrix, warpMatrixType); + } + // return final correlation coefficient + return rho; +} +}; +/* End of file. */ \ No newline at end of file diff --git a/modules/video/test/test_ecc.cpp b/modules/video/test/test_ecc.cpp index 0bb3045dbc..ccbaa65230 100644 --- a/modules/video/test/test_ecc.cpp +++ b/modules/video/test/test_ecc.cpp @@ -45,16 +45,30 @@ namespace opencv_test { namespace { -class CV_ECC_BaseTest : public cvtest::BaseTest { +PARAM_TEST_CASE(Video_ECC, int, bool) +{ + int motionType; + bool usePyramids; + virtual void SetUp() + { + motionType = GET_PARAM(0); + usePyramids = GET_PARAM(1); + } +}; + +class CV_ECC_Test : public cvtest::BaseTest { public: - CV_ECC_BaseTest(); - virtual ~CV_ECC_BaseTest(); + CV_ECC_Test(int motionType, bool usePyramids); + virtual ~CV_ECC_Test(); protected: + int motionType; + double MAX_RMS; // upper bound for RMS error + double computeRMS(const Mat& mat1, const Mat& mat2); bool isMapCorrect(const Mat& mat); - virtual bool test(const Mat) { return true; }; // single test + virtual bool test(const Mat img); bool testAllTypes(const Mat img); // run test for all supported data types (U8, U16, F32, F64) bool testAllChNum(const Mat img); // run test for all supported channels count (gray, RGB) @@ -62,25 +76,28 @@ class CV_ECC_BaseTest : public cvtest::BaseTest { bool checkMap(const Mat& map, const Mat& ground); - double MAX_RMS_ECC; // upper bound for RMS error int ntests; // number of tests per motion type int ECC_iterations; // number of iterations for ECC double ECC_epsilon; // we choose a negative value, so that // ECC_iterations are always executed TermCriteria criteria; + bool usePyramids; // use version of findTransformECC with pyramids }; -CV_ECC_BaseTest::CV_ECC_BaseTest() { - MAX_RMS_ECC = 0.1; - ntests = 3; - ECC_iterations = 50; - ECC_epsilon = -1; //-> negative value means that ECC_Iterations will be executed - criteria = TermCriteria(TermCriteria::COUNT + TermCriteria::EPS, ECC_iterations, ECC_epsilon); -} -CV_ECC_BaseTest::~CV_ECC_BaseTest() {} +CV_ECC_Test::CV_ECC_Test(int a_motionType, bool a_usePyramids) : motionType(a_motionType) + , MAX_RMS(0.1) + , ntests(3) + , ECC_iterations(50) + , ECC_epsilon(-1) + , criteria(TermCriteria::COUNT + TermCriteria::EPS, ECC_iterations, ECC_epsilon) + , usePyramids(a_usePyramids) + {} -bool CV_ECC_BaseTest::isMapCorrect(const Mat& map) { + +CV_ECC_Test::~CV_ECC_Test() {} + +bool CV_ECC_Test::isMapCorrect(const Mat& map) { bool tr = true; float mapVal; for (int i = 0; i < map.rows; i++) @@ -92,7 +109,7 @@ bool CV_ECC_BaseTest::isMapCorrect(const Mat& map) { return tr; } -double CV_ECC_BaseTest::computeRMS(const Mat& mat1, const Mat& mat2) { +double CV_ECC_Test::computeRMS(const Mat& mat1, const Mat& mat2) { CV_Assert(mat1.rows == mat2.rows); CV_Assert(mat1.cols == mat2.cols); @@ -102,13 +119,13 @@ double CV_ECC_BaseTest::computeRMS(const Mat& mat1, const Mat& mat2) { return sqrt(errorMat.dot(errorMat) / (mat1.rows * mat1.cols * mat1.channels())); } -bool CV_ECC_BaseTest::checkMap(const Mat& map, const Mat& ground) { +bool CV_ECC_Test::checkMap(const Mat& map, const Mat& ground) { if (!isMapCorrect(map)) { ts->set_failed_test_info(cvtest::TS::FAIL_INVALID_OUTPUT); return false; } - if (computeRMS(map, ground) > MAX_RMS_ECC) { + if (computeRMS(map, ground) > MAX_RMS) { ts->set_failed_test_info(cvtest::TS::FAIL_BAD_ACCURACY); ts->printf(ts->LOG, "RMS = %f", computeRMS(map, ground)); return false; @@ -116,7 +133,77 @@ bool CV_ECC_BaseTest::checkMap(const Mat& map, const Mat& ground) { return true; } -bool CV_ECC_BaseTest::testAllTypes(const Mat img) { +bool CV_ECC_Test::test(const Mat img) +{ + cv::RNG rng = ts->get_rng(); + + int progress = 0; + + for (int k = 0; k < ntests; k++) { + ts->update_context(this, k, true); + progress = update_progress(progress, k, ntests, 0); + + Mat groundMap; + switch(motionType) + { + case MOTION_TRANSLATION: + groundMap = (Mat_(2, 3) << 1, 0, (rng.uniform(10.f, 20.f)), 0, 1, (rng.uniform(10.f, 20.f))); + break; + case MOTION_EUCLIDEAN: + { + double angle = CV_PI / 30 + CV_PI * rng.uniform((double)-2.f, (double)2.f) / 180; + groundMap = (Mat_(2, 3) << cos(angle), -sin(angle), (rng.uniform(10.f, 20.f)), sin(angle), + cos(angle), (rng.uniform(10.f, 20.f))); + break; + } + case MOTION_AFFINE: + groundMap = (Mat_(2, 3) << (1 - rng.uniform(-0.05f, 0.05f)), (rng.uniform(-0.03f, 0.03f)), + (rng.uniform(10.f, 20.f)), (rng.uniform(-0.03f, 0.03f)), (1 - rng.uniform(-0.05f, 0.05f)), + (rng.uniform(10.f, 20.f))); + break; + case MOTION_HOMOGRAPHY: + groundMap = + (Mat_(3, 3) << (1 - rng.uniform(-0.05f, 0.05f)), (rng.uniform(-0.03f, 0.03f)), + (rng.uniform(10.f, 20.f)), (rng.uniform(-0.03f, 0.03f)), (1 - rng.uniform(-0.05f, 0.05f)), + (rng.uniform(10.f, 20.f)), (rng.uniform(0.0001f, 0.0003f)), (rng.uniform(0.0001f, 0.0003f)), 1.f); + break; + default: + CV_Error(Error::StsBadArg, "Incorrect motion type"); + break; + } + + Mat warpedImage; + + Mat foundMap; + if(motionType == MOTION_HOMOGRAPHY) + { + warpPerspective(img, warpedImage, groundMap, Size(200, 200), INTER_LINEAR + WARP_INVERSE_MAP); + foundMap = Mat::eye(3, 3, CV_32F); + } + else + { + warpAffine(img, warpedImage, groundMap, Size(200, 200), INTER_LINEAR + WARP_INVERSE_MAP); + foundMap = Mat((Mat_(2, 3) << 1, 0, 0, 0, 1, 0)); + } + + + if(usePyramids) + { + ECCParameters params; + params.criteria = criteria; + params.motionType = motionType; + findTransformECCMultiScale(warpedImage, img, foundMap, params); + } + else + findTransformECC(warpedImage, img, foundMap, motionType, criteria); + + if (!checkMap(foundMap, groundMap)) + return false; + } + return true; +} + +bool CV_ECC_Test::testAllTypes(const Mat img) { auto types = {CV_8U, CV_16U, CV_32F, CV_64F}; for (auto type : types) { Mat timg; @@ -127,9 +214,10 @@ bool CV_ECC_BaseTest::testAllTypes(const Mat img) { return true; } -bool CV_ECC_BaseTest::testAllChNum(const Mat img) { - if (!testAllTypes(img)) - return false; +bool CV_ECC_Test::testAllChNum(const Mat img) { + if(!usePyramids) + if (!testAllTypes(img)) + return false; Mat gray; cvtColor(img, gray, COLOR_RGB2GRAY); @@ -139,7 +227,7 @@ bool CV_ECC_BaseTest::testAllChNum(const Mat img) { return true; } -void CV_ECC_BaseTest::run(int) { +void CV_ECC_Test::run(int) { Mat img = imread(string(ts->get_data_path()) + "shared/fruits.png"); if (img.empty()) { ts->printf(ts->LOG, "test image can not be read"); @@ -155,153 +243,22 @@ void CV_ECC_BaseTest::run(int) { ts->set_failed_test_info(cvtest::TS::OK); } -class CV_ECC_Test_Translation : public CV_ECC_BaseTest { - public: - CV_ECC_Test_Translation(); - - protected: - bool test(const Mat); -}; - -CV_ECC_Test_Translation::CV_ECC_Test_Translation() {} - -bool CV_ECC_Test_Translation::test(const Mat testImg) { - cv::RNG rng = ts->get_rng(); - - int progress = 0; - - for (int k = 0; k < ntests; k++) { - ts->update_context(this, k, true); - progress = update_progress(progress, k, ntests, 0); - - Mat translationGround = (Mat_(2, 3) << 1, 0, (rng.uniform(10.f, 20.f)), 0, 1, (rng.uniform(10.f, 20.f))); - - Mat warpedImage; - - warpAffine(testImg, warpedImage, translationGround, Size(200, 200), INTER_LINEAR + WARP_INVERSE_MAP); - - Mat mapTranslation = (Mat_(2, 3) << 1, 0, 0, 0, 1, 0); - - findTransformECC(warpedImage, testImg, mapTranslation, 0, criteria); - - if (!checkMap(mapTranslation, translationGround)) - return false; - } - return true; +TEST_P(Video_ECC, accuracy) { + CV_ECC_Test test(motionType, usePyramids); + test.safe_run(); } -class CV_ECC_Test_Euclidean : public CV_ECC_BaseTest { - public: - CV_ECC_Test_Euclidean(); +INSTANTIATE_TEST_CASE_P(ECCfixtures, Video_ECC, + testing::Values(testing::make_tuple(MOTION_TRANSLATION, false), + testing::make_tuple(MOTION_TRANSLATION, true), + testing::make_tuple(MOTION_EUCLIDEAN, false), + testing::make_tuple(MOTION_EUCLIDEAN, true), + testing::make_tuple(MOTION_AFFINE, false), + testing::make_tuple(MOTION_AFFINE, true), + testing::make_tuple(MOTION_HOMOGRAPHY, false), + testing::make_tuple(MOTION_HOMOGRAPHY, true))); - protected: - bool test(const Mat); -}; - -CV_ECC_Test_Euclidean::CV_ECC_Test_Euclidean() {} - -bool CV_ECC_Test_Euclidean::test(const Mat testImg) { - cv::RNG rng = ts->get_rng(); - - int progress = 0; - for (int k = 0; k < ntests; k++) { - ts->update_context(this, k, true); - progress = update_progress(progress, k, ntests, 0); - - double angle = CV_PI / 30 + CV_PI * rng.uniform((double)-2.f, (double)2.f) / 180; - - Mat euclideanGround = (Mat_(2, 3) << cos(angle), -sin(angle), (rng.uniform(10.f, 20.f)), sin(angle), - cos(angle), (rng.uniform(10.f, 20.f))); - - Mat warpedImage; - - warpAffine(testImg, warpedImage, euclideanGround, Size(200, 200), INTER_LINEAR + WARP_INVERSE_MAP); - - Mat mapEuclidean = (Mat_(2, 3) << 1, 0, 0, 0, 1, 0); - - findTransformECC(warpedImage, testImg, mapEuclidean, 1, criteria); - - if (!checkMap(mapEuclidean, euclideanGround)) - return false; - } - return true; -} - -class CV_ECC_Test_Affine : public CV_ECC_BaseTest { - public: - CV_ECC_Test_Affine(); - - protected: - bool test(const Mat img); -}; - -CV_ECC_Test_Affine::CV_ECC_Test_Affine() {} - -bool CV_ECC_Test_Affine::test(const Mat testImg) { - cv::RNG rng = ts->get_rng(); - - int progress = 0; - for (int k = 0; k < ntests; k++) { - ts->update_context(this, k, true); - progress = update_progress(progress, k, ntests, 0); - - Mat affineGround = (Mat_(2, 3) << (1 - rng.uniform(-0.05f, 0.05f)), (rng.uniform(-0.03f, 0.03f)), - (rng.uniform(10.f, 20.f)), (rng.uniform(-0.03f, 0.03f)), (1 - rng.uniform(-0.05f, 0.05f)), - (rng.uniform(10.f, 20.f))); - - Mat warpedImage; - - warpAffine(testImg, warpedImage, affineGround, Size(200, 200), INTER_LINEAR + WARP_INVERSE_MAP); - - Mat mapAffine = (Mat_(2, 3) << 1, 0, 0, 0, 1, 0); - - findTransformECC(warpedImage, testImg, mapAffine, 2, criteria); - - if (!checkMap(mapAffine, affineGround)) - return false; - } - - return true; -} - -class CV_ECC_Test_Homography : public CV_ECC_BaseTest { - public: - CV_ECC_Test_Homography(); - - protected: - bool test(const Mat testImg); -}; - -CV_ECC_Test_Homography::CV_ECC_Test_Homography() {} - -bool CV_ECC_Test_Homography::test(const Mat testImg) { - cv::RNG rng = ts->get_rng(); - - int progress = 0; - for (int k = 0; k < ntests; k++) { - ts->update_context(this, k, true); - progress = update_progress(progress, k, ntests, 0); - - Mat homoGround = - (Mat_(3, 3) << (1 - rng.uniform(-0.05f, 0.05f)), (rng.uniform(-0.03f, 0.03f)), - (rng.uniform(10.f, 20.f)), (rng.uniform(-0.03f, 0.03f)), (1 - rng.uniform(-0.05f, 0.05f)), - (rng.uniform(10.f, 20.f)), (rng.uniform(0.0001f, 0.0003f)), (rng.uniform(0.0001f, 0.0003f)), 1.f); - - Mat warpedImage; - - warpPerspective(testImg, warpedImage, homoGround, Size(200, 200), INTER_LINEAR + WARP_INVERSE_MAP); - - Mat mapHomography = Mat::eye(3, 3, CV_32F); - - findTransformECC(warpedImage, testImg, mapHomography, 3, criteria); - - if (!checkMap(mapHomography, homoGround)) - return false; - } - return true; -} - -class CV_ECC_Test_Mask : public CV_ECC_BaseTest { +class CV_ECC_Test_Mask : public CV_ECC_Test { public: CV_ECC_Test_Mask(); @@ -309,7 +266,7 @@ class CV_ECC_Test_Mask : public CV_ECC_BaseTest { bool test(const Mat); }; -CV_ECC_Test_Mask::CV_ECC_Test_Mask() {} +CV_ECC_Test_Mask::CV_ECC_Test_Mask():CV_ECC_Test(MOTION_TRANSLATION, false) {} bool CV_ECC_Test_Mask::test(const Mat testImg) { cv::RNG rng = ts->get_rng(); @@ -368,6 +325,58 @@ bool CV_ECC_Test_Mask::test(const Mat testImg) { return true; } +class CV_ECC_BigPictureTest : public CV_ECC_Test { + public: + CV_ECC_BigPictureTest(bool a_maskedVersion) : CV_ECC_Test(MOTION_HOMOGRAPHY, true), maskedVersion(a_maskedVersion) {} + virtual ~CV_ECC_BigPictureTest() {} + protected: + void run(int); + bool maskedVersion; +}; + +void CV_ECC_BigPictureTest::run(int) +{ + Mat largeGray0 = imread(string(ts->get_data_path()) + "shared/halmosh0.jpg", IMREAD_GRAYSCALE); + Mat largeGray1; + Mat roiMask0; + Mat roiMask1; + Mat expectedRes; + bool readError = false; + if(maskedVersion) + { + largeGray1 = imread(string(ts->get_data_path()) + "shared/halmosh2.jpg", IMREAD_GRAYSCALE); + roiMask0 = imread(string(ts->get_data_path()) + "shared/halmosh0mask.png", IMREAD_GRAYSCALE); + roiMask1 = imread(string(ts->get_data_path()) + "shared/halmosh2mask.png", IMREAD_GRAYSCALE); + readError = largeGray0.empty() || largeGray1.empty() || roiMask0.empty() || roiMask1.empty(); + expectedRes = (Mat_(3, 3) << 1.0225, 0.0606, -28.6452, -0.0475, 1.0314, 11.819, 8.21e-06, -3.65e-07, 1); + } + else + { + largeGray1 = imread(string(ts->get_data_path()) + "shared/halmosh1.jpg", IMREAD_GRAYSCALE); + readError = largeGray0.empty() || largeGray1.empty(); + expectedRes = (Mat_(3, 3) << 0.9756, -0.0319, 24.685, 0.013, 0.9808, 7.7453, -2.35e-05, -9.12e-06, 1); + } + + if(readError) + { + ts->printf(ts->LOG, "test image can not be read"); + ts->set_failed_test_info(cvtest::TS::FAIL_INVALID_TEST_DATA); + return; + } + + cv::Mat found = cv::Mat::eye(3, 3, CV_32F); + constexpr int N_ITERS = 20; + constexpr double TERMINATION_EPS = 1e-6; + ECCParameters params; + params.criteria = cv::TermCriteria(cv::TermCriteria::COUNT + cv::TermCriteria::EPS, N_ITERS, TERMINATION_EPS); + params.motionType = MOTION_HOMOGRAPHY; + params.nlevels = 5; + params.itersPerLevel = {5, 10, 300, 300, 1000}; + findTransformECCMultiScale(largeGray0, largeGray1, found, params, roiMask0, roiMask1); + ASSERT_EQ(checkMap(found, expectedRes), true); + ts->set_failed_test_info(cvtest::TS::OK); +} + void testECCProperties(Mat x, float eps) { // The channels are independent Mat y = x.t(); @@ -450,26 +459,17 @@ TEST(Video_ECC_Test_Compute, bug_14657) { EXPECT_NEAR(computeECC(img, img), 1.0f, 1e-5f); } -TEST(Video_ECC_Translation, accuracy) { - CV_ECC_Test_Translation test; - test.safe_run(); -} -TEST(Video_ECC_Euclidean, accuracy) { - CV_ECC_Test_Euclidean test; - test.safe_run(); -} -TEST(Video_ECC_Affine, accuracy) { - CV_ECC_Test_Affine test; - test.safe_run(); -} -TEST(Video_ECC_Homography, accuracy) { - CV_ECC_Test_Homography test; - test.safe_run(); -} TEST(Video_ECC_Mask, accuracy) { CV_ECC_Test_Mask test; test.safe_run(); } - +TEST(Video_ECC_BigMS, accuracy) { + CV_ECC_BigPictureTest test(false); + test.safe_run(); +} +TEST(Video_ECC_BigMS_Mask, accuracy) { + CV_ECC_BigPictureTest test(true); + test.safe_run(); +} } // namespace } // namespace opencv_test From 5d88781f74b49ec6c2f609b273ce05241d8803f1 Mon Sep 17 00:00:00 2001 From: kjg0724 <63202356+kjg0724@users.noreply.github.com> Date: Tue, 21 Apr 2026 18:24:38 +0900 Subject: [PATCH 17/19] Merge pull request #28782 from kjg0724:reduce-simd-optimization MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit core: add platform-specific SIMD for cv::reduce REDUCE_SUM ### Pull Request Rewrite `cv::reduce` SIMD optimization using platform-specific instructions as discussed in #28763. #### Changes **Col reduce (dim=1) — horizontal sum per row:** - ARM DOTPROD: `vdotq_u32()` — single-instruction byte sum, no intermediate flush needed - ARM AArch64 fallback: `vpaddlq` chain (u8→u16→u32) - Intel AVX2: `_mm256_sad_epu8()` — 32 bytes → 4×u64 partial sums per cycle - Intel SSSE3: `_mm_shuffle_epi8` + `_mm_sad_epu8` for cn=4 channel separation - Intel SSE2: `_mm_sad_epu8` for cn=1 - cn=4: hardware deinterleave (`vld4q_u8` on ARM, shuffle+SAD on Intel) **Row reduce (dim=0) — vertical accumulation across rows:** - u16 intermediate accumulator with 256-row flush (halves memory bandwidth vs direct u32) - ARM AArch64: `vaddw_u8` widening add (single instruction vs expand+add) - Intel: unpack + add with u16 buffer **Coverage:** All REDUCE_SUM type combinations — 8U→32S, 8U→32F, 16U→32F, 16S→32F, 32F→32F, 32F→64F, 64F→64F. Non-8U types use universal intrinsics where platform-specific gain is minimal (widening is single-stage, FP has limited alternatives). `REDUCE_AVG` benefits automatically (uses SUM internally). **Dispatch:** `CV_CPU_DISPATCH` with SSE2/AVX2/NEON_DOTPROD/LASX. Fallback hierarchy: DOTPROD → AArch64 NEON → Universal Intrinsics → scalar. #### Benchmark (Apple M3 Pro, MacBook Pro 16-inch 2023) | Path | Type | Speedup vs scalar | |------|------|-------------------| | Col reduce (dim=1) | 8UC1 | **3.0–4.5x** | | Col reduce (dim=1) | 8UC4 | **2.7–5.0x** | | Col reduce (dim=1) | 32FC1 | 1.7–2.3x | | Row reduce (dim=0) | 8UC1 | 1.1–1.5x | Row reduce gains are modest due to memory-bandwidth bound (as expected for vertical accumulation). No regressions on non-target paths (MAX, MIN, SUM2). #### Testing - `opencv_test_core --gtest_filter="*Reduce*:*reduce*"` — 483 tests PASSED - Edge cases: non-aligned dimensions (127×61), cn=2/3 scalar fallback, REDUCE_SUM2 unmodified #### Files changed - `modules/core/src/reduce.simd.hpp` — new, platform-dispatched SIMD implementation - `modules/core/src/reduce.dispatch.cpp` — new, CV_CPU_DISPATCH wrapper - `modules/core/CMakeLists.txt` — add `NEON_DOTPROD` to dispatch list - `modules/core/src/matrix_operations.cpp` — wire dispatch functions into ReduceC/R_Invoker --- modules/core/CMakeLists.txt | 1 + modules/core/src/matrix_operations.cpp | 27 +- modules/core/src/reduce.dispatch.cpp | 30 + modules/core/src/reduce.simd.hpp | 1133 ++++++++++++++++++++++++ 4 files changed, 1188 insertions(+), 3 deletions(-) create mode 100644 modules/core/src/reduce.dispatch.cpp create mode 100644 modules/core/src/reduce.simd.hpp diff --git a/modules/core/CMakeLists.txt b/modules/core/CMakeLists.txt index a341c6a837..0abd534d85 100644 --- a/modules/core/CMakeLists.txt +++ b/modules/core/CMakeLists.txt @@ -12,6 +12,7 @@ ocv_add_dispatched_file(mean SSE2 AVX2 LASX) ocv_add_dispatched_file(merge SSE2 AVX2 LASX) ocv_add_dispatched_file(split SSE2 AVX2 LASX) ocv_add_dispatched_file(sum SSE2 AVX2 LASX) +ocv_add_dispatched_file(reduce SSE2 SSSE3 AVX2 NEON_DOTPROD) ocv_add_dispatched_file(norm SSE2 SSE4_1 AVX AVX2 NEON_DOTPROD LASX) # dispatching for accuracy tests diff --git a/modules/core/src/matrix_operations.cpp b/modules/core/src/matrix_operations.cpp index b99c79a3a2..043820c375 100644 --- a/modules/core/src/matrix_operations.cpp +++ b/modules/core/src/matrix_operations.cpp @@ -341,6 +341,20 @@ cv::Mat cv::Mat::cross(InputArray _m) const namespace cv { +typedef void (*ReduceSumFunc)(const Mat& src, Mat& dst); +ReduceSumFunc getReduceCSumFunc(int sdepth, int ddepth); +ReduceSumFunc getReduceRSumFunc(int sdepth, int ddepth); + +template +struct ReduceR_SIMD +{ + int operator()(const T*, int start, int, WT*, const Op&) const + { + return start; + } +}; + + template class ReduceR_Invoker : public ParallelLoopBody { @@ -364,7 +378,8 @@ public: for( ; --height; ) { src += srcstep; - i = range.start; + ReduceR_SIMD simd_op; + i = simd_op(src, range.start, range.end, buf, op); #if CV_ENABLE_UNROLLED for(; i <= range.end - 4; i += 4 ) { @@ -806,7 +821,10 @@ void cv::reduce(InputArray _src, OutputArray _dst, int dim, int op, int dtype) { if( op == REDUCE_SUM ) { - if(sdepth == CV_8U && ddepth == CV_32S) + ReduceSumFunc simd_func = getReduceRSumFunc(sdepth, ddepth); + if(simd_func) + func = (ReduceFunc)simd_func; + else if(sdepth == CV_8U && ddepth == CV_32S) func = reduceSumR8u32s; else if(sdepth == CV_8U && ddepth == CV_32F) func = reduceSumR8u32f; @@ -881,7 +899,10 @@ void cv::reduce(InputArray _src, OutputArray _dst, int dim, int op, int dtype) { if(op == REDUCE_SUM) { - if(sdepth == CV_8U && ddepth == CV_32S) + ReduceSumFunc simd_func = getReduceCSumFunc(sdepth, ddepth); + if(simd_func) + func = (ReduceFunc)simd_func; + else if(sdepth == CV_8U && ddepth == CV_32S) func = reduceSumC8u32s; else if(sdepth == CV_8U && ddepth == CV_32F) func = reduceSumC8u32f; diff --git a/modules/core/src/reduce.dispatch.cpp b/modules/core/src/reduce.dispatch.cpp new file mode 100644 index 0000000000..0b63560780 --- /dev/null +++ b/modules/core/src/reduce.dispatch.cpp @@ -0,0 +1,30 @@ +// 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 "precomp.hpp" + +#include "reduce.simd.hpp" +#include "reduce.simd_declarations.hpp" + +namespace cv { + +typedef void (*ReduceSumFunc)(const Mat& src, Mat& dst); +ReduceSumFunc getReduceCSumFunc(int sdepth, int ddepth); +ReduceSumFunc getReduceRSumFunc(int sdepth, int ddepth); + +ReduceSumFunc getReduceCSumFunc(int sdepth, int ddepth) +{ + CV_INSTRUMENT_REGION(); + CV_CPU_DISPATCH(getReduceCSumFunc, (sdepth, ddepth), + CV_CPU_DISPATCH_MODES_ALL); +} + +ReduceSumFunc getReduceRSumFunc(int sdepth, int ddepth) +{ + CV_INSTRUMENT_REGION(); + CV_CPU_DISPATCH(getReduceRSumFunc, (sdepth, ddepth), + CV_CPU_DISPATCH_MODES_ALL); +} + +} // namespace cv diff --git a/modules/core/src/reduce.simd.hpp b/modules/core/src/reduce.simd.hpp new file mode 100644 index 0000000000..374853bb29 --- /dev/null +++ b/modules/core/src/reduce.simd.hpp @@ -0,0 +1,1133 @@ +// 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 "precomp.hpp" + +namespace cv { +CV_CPU_OPTIMIZATION_NAMESPACE_BEGIN + +typedef void (*ReduceSumFunc)(const Mat& src, Mat& dst); +ReduceSumFunc getReduceCSumFunc(int sdepth, int ddepth); +ReduceSumFunc getReduceRSumFunc(int sdepth, int ddepth); + +#ifndef CV_CPU_OPTIMIZATION_DECLARATIONS_ONLY + +// ===================================================================== +// Col reduce SUM (dim=1): sum each row into cn output values +// ===================================================================== + +#if (CV_SIMD || CV_SIMD_SCALABLE) + +// --- uchar → int --- +// Uses u16 intermediate accumulator to halve the number of widen operations: +// u8→u16 per iteration, u16→u32 flush every 128 iterations (max u16 value: 128*255=32640 < 65535) +static void reduceColSum_8u32s(const Mat& srcmat, Mat& dstmat) +{ + const int cn = srcmat.channels(); + const int cols = srcmat.cols; + const int width = cols * cn; + + auto body = [&](const Range& range) { + for (int y = range.start; y < range.end; y++) + { + const uchar* src = srcmat.ptr(y); + int* dst = dstmat.ptr(y); + + if (cn == 1) + { +#if defined(CV_NEON_DOT) + // ARMv8.2 DOTPROD: vdotq_u32(acc, src, ones) — single instruction byte sum + // u32 lane max 1020/iter, overflow at ~4.2M iter (67MB) — safe without flush + const uint8x16_t ones = vdupq_n_u8(1); + uint32x4_t v_sum = vdupq_n_u32(0); + int x = 0; + for (; x <= width - 16; x += 16) + v_sum = vdotq_u32(v_sum, vld1q_u8(src + x), ones); + int total = (int)vaddvq_u32(v_sum); + for (; x < width; x++) + total += (int)src[x]; + dst[0] = total; +#elif CV_NEON && (defined(__aarch64__) || defined(_M_ARM64)) + // AArch64 fallback: vpaddlq chain + uint32x4_t v_sum = vdupq_n_u32(0); + int x = 0; + for (; x <= width - 16; x += 16) + v_sum = vaddq_u32(v_sum, vpaddlq_u16(vpaddlq_u8(vld1q_u8(src + x)))); + int total = (int)vaddvq_u32(v_sum); + for (; x < width; x++) + total += (int)src[x]; + dst[0] = total; +#elif CV_AVX2 + // Intel AVX2: _mm256_sad_epu8 — 32 bytes → 4×u64 sum, 1 cycle + const __m256i zero = _mm256_setzero_si256(); + __m256i vsum = zero; + int x = 0; + for (; x <= width - 32; x += 32) + vsum = _mm256_add_epi64(vsum, _mm256_sad_epu8( + _mm256_loadu_si256((const __m256i*)(src + x)), zero)); + // 8-byte tail: reduce max scalar tail from 31 to 7 + // Use _mm256_inserti128_si256 instead of _mm256_castsi128_si256: + // the latter leaves upper 128 bits undefined per Intel spec, + // causing MSVC to expose stale register data to _mm256_sad_epu8. + for (; x <= width - 8; x += 8) + vsum = _mm256_add_epi64(vsum, _mm256_sad_epu8( + _mm256_inserti128_si256(_mm256_setzero_si256(), _mm_loadl_epi64((const __m128i*)(src + x)), 0), zero)); + __m128i lo128 = _mm256_castsi256_si128(vsum); + __m128i hi128 = _mm256_extracti128_si256(vsum, 1); + __m128i s = _mm_add_epi64(lo128, hi128); + s = _mm_add_epi64(s, _mm_unpackhi_epi64(s, s)); + int total = _mm_cvtsi128_si32(s); + for (; x < width; x++) + total += (int)src[x]; + dst[0] = total; +#elif CV_SSE2 + // Intel SSE2: _mm_sad_epu8 — 16 bytes → 2×u64 sum + const __m128i zero = _mm_setzero_si128(); + __m128i vsum = zero; + int x = 0; + for (; x <= width - 16; x += 16) + vsum = _mm_add_epi64(vsum, _mm_sad_epu8( + _mm_loadu_si128((const __m128i*)(src + x)), zero)); + __m128i s = _mm_add_epi64(vsum, _mm_unpackhi_epi64(vsum, vsum)); + int total = _mm_cvtsi128_si32(s); + for (; x < width; x++) + total += (int)src[x]; + dst[0] = total; +#else + const int vlanes8 = VTraits::vlanes(); + v_uint32 v_sum = vx_setzero_u32(); + v_uint16 v_sum16 = vx_setzero_u16(); + int x = 0, batch = 0; + const int flush_at = 128; + + for (; x <= width - vlanes8; x += vlanes8) + { + v_uint16 lo, hi; + v_expand(vx_load(src + x), lo, hi); + v_sum16 = v_add(v_sum16, v_add(lo, hi)); + + if (++batch >= flush_at) + { + v_uint32 a, b; + v_expand(v_sum16, a, b); + v_sum = v_add(v_sum, v_add(a, b)); + v_sum16 = vx_setzero_u16(); + batch = 0; + } + } + { + v_uint32 a, b; + v_expand(v_sum16, a, b); + v_sum = v_add(v_sum, v_add(a, b)); + } + + int total = (int)v_reduce_sum(v_sum); + for (; x < width; x++) + total += (int)src[x]; + dst[0] = total; +#endif + } + else if (cn == 4) + { +#if defined(CV_NEON_DOT) + // ARMv8.2 DOTPROD: vld4q_u8 deinterleave + vdotq_u32 per-channel + // vld4q_u8 reads 16 pixels (64 bytes), deinterleaves into 4 × u8x16 + // vdotq_u32 sums each channel's 16 bytes into 4 partial u32 sums + const uint8x16_t ones = vdupq_n_u8(1); + uint32x4_t acc0 = vdupq_n_u32(0), acc1 = vdupq_n_u32(0); + uint32x4_t acc2 = vdupq_n_u32(0), acc3 = vdupq_n_u32(0); + int x = 0; + for (; x <= cols - 16; x += 16) + { + uint8x16x4_t v = vld4q_u8(src + x * 4); + acc0 = vdotq_u32(acc0, v.val[0], ones); + acc1 = vdotq_u32(acc1, v.val[1], ones); + acc2 = vdotq_u32(acc2, v.val[2], ones); + acc3 = vdotq_u32(acc3, v.val[3], ones); + } + dst[0] = (int)vaddvq_u32(acc0); + dst[1] = (int)vaddvq_u32(acc1); + dst[2] = (int)vaddvq_u32(acc2); + dst[3] = (int)vaddvq_u32(acc3); + for (; x < cols; x++) + { + dst[0] += (int)src[x * 4]; + dst[1] += (int)src[x * 4 + 1]; + dst[2] += (int)src[x * 4 + 2]; + dst[3] += (int)src[x * 4 + 3]; + } +#elif CV_NEON && (defined(__aarch64__) || defined(_M_ARM64)) + // AArch64 fallback: vld4q_u8 deinterleave + vaddlvq_u8 per-channel + int sums[4] = {0, 0, 0, 0}; + int x = 0; + for (; x <= cols - 16; x += 16) + { + uint8x16x4_t v = vld4q_u8(src + x * 4); + sums[0] += (int)vaddlvq_u8(v.val[0]); + sums[1] += (int)vaddlvq_u8(v.val[1]); + sums[2] += (int)vaddlvq_u8(v.val[2]); + sums[3] += (int)vaddlvq_u8(v.val[3]); + } + dst[0] = sums[0]; dst[1] = sums[1]; + dst[2] = sums[2]; dst[3] = sums[3]; + for (; x < cols; x++) + { + dst[0] += (int)src[x * 4]; + dst[1] += (int)src[x * 4 + 1]; + dst[2] += (int)src[x * 4 + 2]; + dst[3] += (int)src[x * 4 + 3]; + } +#elif CV_AVX2 + // Intel AVX2: shuffle + unpack to group channels, full-register SAD + const __m256i zero = _mm256_setzero_si256(); + const __m256i shuf_mask = _mm256_setr_epi8( + 0,4,8,12, 1,5,9,13, 2,6,10,14, 3,7,11,15, + 0,4,8,12, 1,5,9,13, 2,6,10,14, 3,7,11,15); + int x = 0; + // 16 pixels/iter with unpack (full-register SAD) + { + __m256i sum_bg = zero, sum_ra = zero; + for (; x <= cols - 16; x += 16) + { + __m256i v0 = _mm256_loadu_si256((const __m256i*)(src + x * 4)); + __m256i v1 = _mm256_loadu_si256((const __m256i*)(src + x * 4 + 32)); + __m256i g0 = _mm256_shuffle_epi8(v0, shuf_mask); + __m256i g1 = _mm256_shuffle_epi8(v1, shuf_mask); + __m256i bg = _mm256_unpacklo_epi32(g0, g1); + __m256i ra = _mm256_unpackhi_epi32(g0, g1); + sum_bg = _mm256_add_epi64(sum_bg, _mm256_sad_epu8(bg, zero)); + sum_ra = _mm256_add_epi64(sum_ra, _mm256_sad_epu8(ra, zero)); + } + __m128i bg_lo = _mm256_castsi256_si128(sum_bg); + __m128i bg_hi = _mm256_extracti128_si256(sum_bg, 1); + __m128i bg_s = _mm_add_epi64(bg_lo, bg_hi); + __m128i ra_lo = _mm256_castsi256_si128(sum_ra); + __m128i ra_hi = _mm256_extracti128_si256(sum_ra, 1); + __m128i ra_s = _mm_add_epi64(ra_lo, ra_hi); + dst[0] = _mm_cvtsi128_si32(bg_s); // B + dst[1] = _mm_cvtsi128_si32(_mm_srli_si128(bg_s, 8)); // G + dst[2] = _mm_cvtsi128_si32(ra_s); // R + dst[3] = _mm_cvtsi128_si32(_mm_srli_si128(ra_s, 8)); // A + } + // 8-pixel tail with and/srli + { + const __m256i mask_lo32 = _mm256_set1_epi64x(0x00000000FFFFFFFFLL); + __m256i sum_br = zero, sum_ga = zero; + for (; x <= cols - 8; x += 8) + { + __m256i v = _mm256_loadu_si256((const __m256i*)(src + x * 4)); + __m256i grouped = _mm256_shuffle_epi8(v, shuf_mask); + __m256i br = _mm256_and_si256(grouped, mask_lo32); + __m256i ga = _mm256_srli_epi64(grouped, 32); + sum_br = _mm256_add_epi64(sum_br, _mm256_sad_epu8(br, zero)); + sum_ga = _mm256_add_epi64(sum_ga, _mm256_sad_epu8(ga, zero)); + } + __m128i br_lo = _mm256_castsi256_si128(sum_br); + __m128i br_hi = _mm256_extracti128_si256(sum_br, 1); + __m128i br_s = _mm_add_epi64(br_lo, br_hi); + __m128i ga_lo = _mm256_castsi256_si128(sum_ga); + __m128i ga_hi = _mm256_extracti128_si256(sum_ga, 1); + __m128i ga_s = _mm_add_epi64(ga_lo, ga_hi); + dst[0] += _mm_cvtsi128_si32(br_s); + dst[1] += _mm_cvtsi128_si32(ga_s); + dst[2] += _mm_cvtsi128_si32(_mm_srli_si128(br_s, 8)); + dst[3] += _mm_cvtsi128_si32(_mm_srli_si128(ga_s, 8)); + } + for (; x < cols; x++) + { + dst[0] += (int)src[x * 4]; + dst[1] += (int)src[x * 4 + 1]; + dst[2] += (int)src[x * 4 + 2]; + dst[3] += (int)src[x * 4 + 3]; + } +#elif CV_SSSE3 + // Intel SSSE3: _mm_shuffle_epi8 deinterleave + _mm_sad_epu8 + // 16 bytes = 4 pixels × 4 channels per iteration + const __m128i zero = _mm_setzero_si128(); + const __m128i shuf_mask = _mm_setr_epi8( + 0,4,8,12, 1,5,9,13, 2,6,10,14, 3,7,11,15); + const __m128i mask_lo32 = _mm_set1_epi64x(0x00000000FFFFFFFFLL); + __m128i sum_br = zero, sum_ga = zero; + int x = 0; + for (; x <= cols - 4; x += 4) + { + __m128i v = _mm_loadu_si128((const __m128i*)(src + x * 4)); + __m128i grouped = _mm_shuffle_epi8(v, shuf_mask); + __m128i br = _mm_and_si128(grouped, mask_lo32); + __m128i ga = _mm_srli_epi64(grouped, 32); + sum_br = _mm_add_epi64(sum_br, _mm_sad_epu8(br, zero)); + sum_ga = _mm_add_epi64(sum_ga, _mm_sad_epu8(ga, zero)); + } + dst[0] = _mm_cvtsi128_si32(sum_br); + dst[1] = _mm_cvtsi128_si32(sum_ga); + dst[2] = _mm_cvtsi128_si32(_mm_srli_si128(sum_br, 8)); + dst[3] = _mm_cvtsi128_si32(_mm_srli_si128(sum_ga, 8)); + for (; x < cols; x++) + { + dst[0] += (int)src[x * 4]; + dst[1] += (int)src[x * 4 + 1]; + dst[2] += (int)src[x * 4 + 2]; + dst[3] += (int)src[x * 4 + 3]; + } +#else + const int vlanes8 = VTraits::vlanes(); + v_uint32 vsum0 = vx_setzero_u32(), vsum1 = vx_setzero_u32(); + v_uint32 vsum2 = vx_setzero_u32(), vsum3 = vx_setzero_u32(); + v_uint16 vs16_0 = vx_setzero_u16(), vs16_1 = vx_setzero_u16(); + v_uint16 vs16_2 = vx_setzero_u16(), vs16_3 = vx_setzero_u16(); + int x = 0, batch = 0; + const int flush_at = 128; + + for (; x <= cols - vlanes8; x += vlanes8) + { + v_uint8 ch0, ch1, ch2, ch3; + v_load_deinterleave(src + x * 4, ch0, ch1, ch2, ch3); + + v_uint16 lo, hi; + v_expand(ch0, lo, hi); vs16_0 = v_add(vs16_0, v_add(lo, hi)); + v_expand(ch1, lo, hi); vs16_1 = v_add(vs16_1, v_add(lo, hi)); + v_expand(ch2, lo, hi); vs16_2 = v_add(vs16_2, v_add(lo, hi)); + v_expand(ch3, lo, hi); vs16_3 = v_add(vs16_3, v_add(lo, hi)); + + if (++batch >= flush_at) + { + v_uint32 a, b; + v_expand(vs16_0, a, b); vsum0 = v_add(vsum0, v_add(a, b)); + v_expand(vs16_1, a, b); vsum1 = v_add(vsum1, v_add(a, b)); + v_expand(vs16_2, a, b); vsum2 = v_add(vsum2, v_add(a, b)); + v_expand(vs16_3, a, b); vsum3 = v_add(vsum3, v_add(a, b)); + vs16_0 = vx_setzero_u16(); vs16_1 = vx_setzero_u16(); + vs16_2 = vx_setzero_u16(); vs16_3 = vx_setzero_u16(); + batch = 0; + } + } + { + v_uint32 a, b; + v_expand(vs16_0, a, b); vsum0 = v_add(vsum0, v_add(a, b)); + v_expand(vs16_1, a, b); vsum1 = v_add(vsum1, v_add(a, b)); + v_expand(vs16_2, a, b); vsum2 = v_add(vsum2, v_add(a, b)); + v_expand(vs16_3, a, b); vsum3 = v_add(vsum3, v_add(a, b)); + } + + dst[0] = (int)v_reduce_sum(vsum0); + dst[1] = (int)v_reduce_sum(vsum1); + dst[2] = (int)v_reduce_sum(vsum2); + dst[3] = (int)v_reduce_sum(vsum3); + + for (; x < cols; x++) + { + dst[0] += (int)src[x * 4]; + dst[1] += (int)src[x * 4 + 1]; + dst[2] += (int)src[x * 4 + 2]; + dst[3] += (int)src[x * 4 + 3]; + } +#endif + } + else + { + // generic cn: scalar fallback + for (int c = 0; c < cn; c++) + dst[c] = (int)src[c]; + for (int x = cn; x < width; x += cn) + for (int c = 0; c < cn; c++) + dst[c] += (int)src[x + c]; + } + } + }; + + parallel_for_(Range(0, srcmat.rows), body); + v_cleanup(); +} + +// --- uchar → float --- +static void reduceColSum_8u32f(const Mat& srcmat, Mat& dstmat) +{ + const int cn = srcmat.channels(); + const int cols = srcmat.cols; + const int width = cols * cn; + + auto body = [&](const Range& range) { + for (int y = range.start; y < range.end; y++) + { + const uchar* src = srcmat.ptr(y); + float* dst = dstmat.ptr(y); + + // compute in int, then convert to float + AutoBuffer ibuf(cn); + int* sums = ibuf.data(); + for (int c = 0; c < cn; c++) + sums[c] = 0; + + if (cn == 1) + { +#if defined(CV_NEON_DOT) + const uint8x16_t ones = vdupq_n_u8(1); + uint32x4_t v_sum = vdupq_n_u32(0); + int x = 0; + for (; x <= width - 16; x += 16) + v_sum = vdotq_u32(v_sum, vld1q_u8(src + x), ones); + sums[0] = (int)vaddvq_u32(v_sum); + for (; x < width; x++) + sums[0] += (int)src[x]; +#elif CV_NEON && (defined(__aarch64__) || defined(_M_ARM64)) + uint32x4_t v_sum = vdupq_n_u32(0); + int x = 0; + for (; x <= width - 16; x += 16) + v_sum = vaddq_u32(v_sum, vpaddlq_u16(vpaddlq_u8(vld1q_u8(src + x)))); + sums[0] = (int)vaddvq_u32(v_sum); + for (; x < width; x++) + sums[0] += (int)src[x]; +#elif CV_AVX2 + const __m256i zero = _mm256_setzero_si256(); + __m256i vsum = zero; + int x = 0; + for (; x <= width - 32; x += 32) + vsum = _mm256_add_epi64(vsum, _mm256_sad_epu8( + _mm256_loadu_si256((const __m256i*)(src + x)), zero)); + // 8-byte tail: reduce max scalar tail from 31 to 7 + // Use _mm256_inserti128_si256 instead of _mm256_castsi128_si256: + // the latter leaves upper 128 bits undefined per Intel spec, + // causing MSVC to expose stale register data to _mm256_sad_epu8. + for (; x <= width - 8; x += 8) + vsum = _mm256_add_epi64(vsum, _mm256_sad_epu8( + _mm256_inserti128_si256(_mm256_setzero_si256(), _mm_loadl_epi64((const __m128i*)(src + x)), 0), zero)); + __m128i lo128 = _mm256_castsi256_si128(vsum); + __m128i hi128 = _mm256_extracti128_si256(vsum, 1); + __m128i s = _mm_add_epi64(lo128, hi128); + s = _mm_add_epi64(s, _mm_unpackhi_epi64(s, s)); + sums[0] = _mm_cvtsi128_si32(s); + for (; x < width; x++) + sums[0] += (int)src[x]; +#elif CV_SSE2 + const __m128i zero = _mm_setzero_si128(); + __m128i vsum = zero; + int x = 0; + for (; x <= width - 16; x += 16) + vsum = _mm_add_epi64(vsum, _mm_sad_epu8( + _mm_loadu_si128((const __m128i*)(src + x)), zero)); + __m128i s = _mm_add_epi64(vsum, _mm_unpackhi_epi64(vsum, vsum)); + sums[0] = _mm_cvtsi128_si32(s); + for (; x < width; x++) + sums[0] += (int)src[x]; +#else + const int vlanes8 = VTraits::vlanes(); + v_uint32 v_sum = vx_setzero_u32(); + v_uint16 v_sum16 = vx_setzero_u16(); + int x = 0, batch = 0; + + for (; x <= width - vlanes8; x += vlanes8) + { + v_uint16 lo, hi; + v_expand(vx_load(src + x), lo, hi); + v_sum16 = v_add(v_sum16, v_add(lo, hi)); + if (++batch >= 128) + { + v_uint32 a, b; + v_expand(v_sum16, a, b); + v_sum = v_add(v_sum, v_add(a, b)); + v_sum16 = vx_setzero_u16(); + batch = 0; + } + } + v_uint32 a, b; + v_expand(v_sum16, a, b); + v_sum = v_add(v_sum, v_add(a, b)); + + sums[0] = (int)v_reduce_sum(v_sum); + for (; x < width; x++) + sums[0] += (int)src[x]; +#endif + } + else if (cn == 4) + { +#if defined(CV_NEON_DOT) + const uint8x16_t ones = vdupq_n_u8(1); + uint32x4_t acc0 = vdupq_n_u32(0), acc1 = vdupq_n_u32(0); + uint32x4_t acc2 = vdupq_n_u32(0), acc3 = vdupq_n_u32(0); + int x = 0; + for (; x <= cols - 16; x += 16) + { + uint8x16x4_t v = vld4q_u8(src + x * 4); + acc0 = vdotq_u32(acc0, v.val[0], ones); + acc1 = vdotq_u32(acc1, v.val[1], ones); + acc2 = vdotq_u32(acc2, v.val[2], ones); + acc3 = vdotq_u32(acc3, v.val[3], ones); + } + sums[0] = (int)vaddvq_u32(acc0); + sums[1] = (int)vaddvq_u32(acc1); + sums[2] = (int)vaddvq_u32(acc2); + sums[3] = (int)vaddvq_u32(acc3); + for (; x < cols; x++) + { + sums[0] += (int)src[x * 4]; + sums[1] += (int)src[x * 4 + 1]; + sums[2] += (int)src[x * 4 + 2]; + sums[3] += (int)src[x * 4 + 3]; + } +#elif CV_NEON && (defined(__aarch64__) || defined(_M_ARM64)) + int x = 0; + for (; x <= cols - 16; x += 16) + { + uint8x16x4_t v = vld4q_u8(src + x * 4); + sums[0] += (int)vaddlvq_u8(v.val[0]); + sums[1] += (int)vaddlvq_u8(v.val[1]); + sums[2] += (int)vaddlvq_u8(v.val[2]); + sums[3] += (int)vaddlvq_u8(v.val[3]); + } + for (; x < cols; x++) + { + sums[0] += (int)src[x * 4]; + sums[1] += (int)src[x * 4 + 1]; + sums[2] += (int)src[x * 4 + 2]; + sums[3] += (int)src[x * 4 + 3]; + } +#elif CV_AVX2 + // Intel AVX2: shuffle + unpack to group channels, full-register SAD + const __m256i zero = _mm256_setzero_si256(); + const __m256i shuf_mask = _mm256_setr_epi8( + 0,4,8,12, 1,5,9,13, 2,6,10,14, 3,7,11,15, + 0,4,8,12, 1,5,9,13, 2,6,10,14, 3,7,11,15); + int x = 0; + // 16 pixels/iter with unpack (full-register SAD) + { + __m256i sum_bg = zero, sum_ra = zero; + for (; x <= cols - 16; x += 16) + { + __m256i v0 = _mm256_loadu_si256((const __m256i*)(src + x * 4)); + __m256i v1 = _mm256_loadu_si256((const __m256i*)(src + x * 4 + 32)); + __m256i g0 = _mm256_shuffle_epi8(v0, shuf_mask); + __m256i g1 = _mm256_shuffle_epi8(v1, shuf_mask); + __m256i bg = _mm256_unpacklo_epi32(g0, g1); + __m256i ra = _mm256_unpackhi_epi32(g0, g1); + sum_bg = _mm256_add_epi64(sum_bg, _mm256_sad_epu8(bg, zero)); + sum_ra = _mm256_add_epi64(sum_ra, _mm256_sad_epu8(ra, zero)); + } + __m128i bg_lo = _mm256_castsi256_si128(sum_bg); + __m128i bg_hi = _mm256_extracti128_si256(sum_bg, 1); + __m128i bg_s = _mm_add_epi64(bg_lo, bg_hi); + __m128i ra_lo = _mm256_castsi256_si128(sum_ra); + __m128i ra_hi = _mm256_extracti128_si256(sum_ra, 1); + __m128i ra_s = _mm_add_epi64(ra_lo, ra_hi); + sums[0] = _mm_cvtsi128_si32(bg_s); // B + sums[1] = _mm_cvtsi128_si32(_mm_srli_si128(bg_s, 8)); // G + sums[2] = _mm_cvtsi128_si32(ra_s); // R + sums[3] = _mm_cvtsi128_si32(_mm_srli_si128(ra_s, 8)); // A + } + // 8-pixel tail with and/srli + { + const __m256i mask_lo32 = _mm256_set1_epi64x(0x00000000FFFFFFFFLL); + __m256i sum_br = zero, sum_ga = zero; + for (; x <= cols - 8; x += 8) + { + __m256i v = _mm256_loadu_si256((const __m256i*)(src + x * 4)); + __m256i grouped = _mm256_shuffle_epi8(v, shuf_mask); + __m256i br = _mm256_and_si256(grouped, mask_lo32); + __m256i ga = _mm256_srli_epi64(grouped, 32); + sum_br = _mm256_add_epi64(sum_br, _mm256_sad_epu8(br, zero)); + sum_ga = _mm256_add_epi64(sum_ga, _mm256_sad_epu8(ga, zero)); + } + __m128i br_lo = _mm256_castsi256_si128(sum_br); + __m128i br_hi = _mm256_extracti128_si256(sum_br, 1); + __m128i br_s = _mm_add_epi64(br_lo, br_hi); + __m128i ga_lo = _mm256_castsi256_si128(sum_ga); + __m128i ga_hi = _mm256_extracti128_si256(sum_ga, 1); + __m128i ga_s = _mm_add_epi64(ga_lo, ga_hi); + sums[0] += _mm_cvtsi128_si32(br_s); + sums[1] += _mm_cvtsi128_si32(ga_s); + sums[2] += _mm_cvtsi128_si32(_mm_srli_si128(br_s, 8)); + sums[3] += _mm_cvtsi128_si32(_mm_srli_si128(ga_s, 8)); + } + for (; x < cols; x++) + { + sums[0] += (int)src[x * 4]; + sums[1] += (int)src[x * 4 + 1]; + sums[2] += (int)src[x * 4 + 2]; + sums[3] += (int)src[x * 4 + 3]; + } +#elif CV_SSSE3 + // Intel SSSE3: _mm_shuffle_epi8 deinterleave + _mm_sad_epu8 + const __m128i zero = _mm_setzero_si128(); + const __m128i shuf_mask = _mm_setr_epi8( + 0,4,8,12, 1,5,9,13, 2,6,10,14, 3,7,11,15); + const __m128i mask_lo32 = _mm_set1_epi64x(0x00000000FFFFFFFFLL); + __m128i sum_br = zero, sum_ga = zero; + int x = 0; + for (; x <= cols - 4; x += 4) + { + __m128i v = _mm_loadu_si128((const __m128i*)(src + x * 4)); + __m128i grouped = _mm_shuffle_epi8(v, shuf_mask); + __m128i br = _mm_and_si128(grouped, mask_lo32); + __m128i ga = _mm_srli_epi64(grouped, 32); + sum_br = _mm_add_epi64(sum_br, _mm_sad_epu8(br, zero)); + sum_ga = _mm_add_epi64(sum_ga, _mm_sad_epu8(ga, zero)); + } + sums[0] = _mm_cvtsi128_si32(sum_br); + sums[1] = _mm_cvtsi128_si32(sum_ga); + sums[2] = _mm_cvtsi128_si32(_mm_srli_si128(sum_br, 8)); + sums[3] = _mm_cvtsi128_si32(_mm_srli_si128(sum_ga, 8)); + for (; x < cols; x++) + { + sums[0] += (int)src[x * 4]; + sums[1] += (int)src[x * 4 + 1]; + sums[2] += (int)src[x * 4 + 2]; + sums[3] += (int)src[x * 4 + 3]; + } +#else + const int vlanes8 = VTraits::vlanes(); + v_uint32 vsum0 = vx_setzero_u32(), vsum1 = vx_setzero_u32(); + v_uint32 vsum2 = vx_setzero_u32(), vsum3 = vx_setzero_u32(); + v_uint16 vs16_0 = vx_setzero_u16(), vs16_1 = vx_setzero_u16(); + v_uint16 vs16_2 = vx_setzero_u16(), vs16_3 = vx_setzero_u16(); + int x = 0, batch = 0; + const int flush_at = 128; + + for (; x <= cols - vlanes8; x += vlanes8) + { + v_uint8 ch0, ch1, ch2, ch3; + v_load_deinterleave(src + x * 4, ch0, ch1, ch2, ch3); + + v_uint16 lo, hi; + v_expand(ch0, lo, hi); vs16_0 = v_add(vs16_0, v_add(lo, hi)); + v_expand(ch1, lo, hi); vs16_1 = v_add(vs16_1, v_add(lo, hi)); + v_expand(ch2, lo, hi); vs16_2 = v_add(vs16_2, v_add(lo, hi)); + v_expand(ch3, lo, hi); vs16_3 = v_add(vs16_3, v_add(lo, hi)); + + if (++batch >= flush_at) + { + v_uint32 a, b; + v_expand(vs16_0, a, b); vsum0 = v_add(vsum0, v_add(a, b)); + v_expand(vs16_1, a, b); vsum1 = v_add(vsum1, v_add(a, b)); + v_expand(vs16_2, a, b); vsum2 = v_add(vsum2, v_add(a, b)); + v_expand(vs16_3, a, b); vsum3 = v_add(vsum3, v_add(a, b)); + vs16_0 = vx_setzero_u16(); vs16_1 = vx_setzero_u16(); + vs16_2 = vx_setzero_u16(); vs16_3 = vx_setzero_u16(); + batch = 0; + } + } + { + v_uint32 a, b; + v_expand(vs16_0, a, b); vsum0 = v_add(vsum0, v_add(a, b)); + v_expand(vs16_1, a, b); vsum1 = v_add(vsum1, v_add(a, b)); + v_expand(vs16_2, a, b); vsum2 = v_add(vsum2, v_add(a, b)); + v_expand(vs16_3, a, b); vsum3 = v_add(vsum3, v_add(a, b)); + } + + sums[0] = (int)v_reduce_sum(vsum0); + sums[1] = (int)v_reduce_sum(vsum1); + sums[2] = (int)v_reduce_sum(vsum2); + sums[3] = (int)v_reduce_sum(vsum3); + + for (; x < cols; x++) + { + sums[0] += (int)src[x * 4]; + sums[1] += (int)src[x * 4 + 1]; + sums[2] += (int)src[x * 4 + 2]; + sums[3] += (int)src[x * 4 + 3]; + } +#endif + } + else + { + for (int x = 0, c = 0; x < width; x++, c = (c+1)&-(c < cn-1)) + sums[c] += (int)src[x]; + } + for (int c = 0; c < cn; c++) + dst[c] = (float)sums[c]; + } + }; + parallel_for_(Range(0, srcmat.rows), body); + v_cleanup(); +} + +// --- 16-bit (ushort/short) → float --- +template +static void reduceColSum_16_32f(const Mat& srcmat, Mat& dstmat) +{ + const int cn = srcmat.channels(); + const int width = srcmat.cols * cn; + const int vlanes = VTraits::vlanes(); + + auto body = [&](const Range& range) { + for (int y = range.start; y < range.end; y++) + { + const SrcT* src = srcmat.ptr(y); + float* dst = dstmat.ptr(y); + + if (cn == 1) + { + VecDst v_sum0 = v_setzero_(), v_sum1 = v_setzero_(); + int x = 0; + for (; x <= width - vlanes; x += vlanes) + { + VecDst lo, hi; + v_expand(vx_load(src + x), lo, hi); + v_sum0 = v_add(v_sum0, lo); + v_sum1 = v_add(v_sum1, hi); + } + dst[0] = (float)(v_reduce_sum(v_sum0) + v_reduce_sum(v_sum1)); + for (; x < width; x++) + dst[0] += (float)src[x]; + } + else + { + for (int c = 0; c < cn; c++) + dst[c] = 0; + for (int x = 0, c = 0; x < width; x++, c = (c+1)&-(c < cn-1)) + dst[c] += (float)src[x]; + } + } + }; + parallel_for_(Range(0, srcmat.rows), body); + v_cleanup(); +} + +static void reduceColSum_16u32f(const Mat& srcmat, Mat& dstmat) +{ reduceColSum_16_32f(srcmat, dstmat); } + +static void reduceColSum_16s32f(const Mat& srcmat, Mat& dstmat) +{ reduceColSum_16_32f(srcmat, dstmat); } + +// --- float/double → same type --- +template +static void reduceColSum_FP(const Mat& srcmat, Mat& dstmat) +{ + const int cn = srcmat.channels(); + const int width = srcmat.cols * cn; + const int vlanes = VTraits::vlanes(); + + auto body = [&](const Range& range) { + for (int y = range.start; y < range.end; y++) + { + const T* src = srcmat.ptr(y); + T* dst = dstmat.ptr(y); + + if (cn == 1) + { + VecT s0 = v_setzero_(), s1 = v_setzero_(); + VecT s2 = v_setzero_(), s3 = v_setzero_(); + int x = 0; + for (; x <= width - vlanes * 4; x += vlanes * 4) + { + s0 = v_add(s0, vx_load(src + x)); + s1 = v_add(s1, vx_load(src + x + vlanes)); + s2 = v_add(s2, vx_load(src + x + vlanes * 2)); + s3 = v_add(s3, vx_load(src + x + vlanes * 3)); + } + s0 = v_add(v_add(s0, s1), v_add(s2, s3)); + for (; x <= width - vlanes; x += vlanes) + s0 = v_add(s0, vx_load(src + x)); + T total = (T)v_reduce_sum(s0); + for (; x < width; x++) + total += src[x]; + dst[0] = total; + } + else + { + for (int c = 0; c < cn; c++) + dst[c] = 0; + for (int x = 0, c = 0; x < width; x++, c = (c+1)&-(c < cn-1)) + dst[c] += src[x]; + } + } + }; + parallel_for_(Range(0, srcmat.rows), body); + v_cleanup(); +} + +static void reduceColSum_32f32f(const Mat& srcmat, Mat& dstmat) +{ reduceColSum_FP(srcmat, dstmat); } + +#if CV_SIMD_64F +// --- float → double --- +static void reduceColSum_32f64f(const Mat& srcmat, Mat& dstmat) +{ + const int cn = srcmat.channels(); + const int width = srcmat.cols * cn; + const int vlanes32 = VTraits::vlanes(); + + auto body = [&](const Range& range) { + for (int y = range.start; y < range.end; y++) + { + const float* src = srcmat.ptr(y); + double* dst = dstmat.ptr(y); + + if (cn == 1) + { + v_float64 v_sum0 = vx_setzero_f64(); + v_float64 v_sum1 = vx_setzero_f64(); + int x = 0; + for (; x <= width - vlanes32; x += vlanes32) + { + v_float32 v = vx_load(src + x); + v_sum0 = v_add(v_sum0, v_cvt_f64(v)); + v_sum1 = v_add(v_sum1, v_cvt_f64_high(v)); + } + v_float64 v_total = v_add(v_sum0, v_sum1); + + double total = v_reduce_sum(v_total); + for (; x < width; x++) + total += (double)src[x]; + dst[0] = total; + } + else + { + for (int c = 0; c < cn; c++) + dst[c] = 0; + for (int x = 0, c = 0; x < width; x++, c = (c+1)&-(c < cn-1)) + dst[c] += (double)src[x]; + } + } + }; + parallel_for_(Range(0, srcmat.rows), body); + v_cleanup(); +} + +static void reduceColSum_64f64f(const Mat& srcmat, Mat& dstmat) +{ reduceColSum_FP(srcmat, dstmat); } +#endif // CV_SIMD_64F + +// ===================================================================== +// Row reduce SUM (dim=0): sum each column across all rows +// ===================================================================== + +// --- uchar → int --- +// Uses u16 intermediate accumulator for vertical sum: accumulate u8→u16 per row, +// flush u16→u32 every 256 rows (max u16 value: 256*255=65280 < 65535) +static void reduceRowSum_8u32s(const Mat& srcmat, Mat& dstmat) +{ + const int width_cn = srcmat.cols * srcmat.channels(); + const int height = srcmat.rows; + const int vlanes16 = VTraits::vlanes(); + const int vlanes32 = VTraits::vlanes(); + + auto body = [&](const Range& range) { + const int start = range.start; + const int end = range.end; + const int len = end - start; + + AutoBuffer buf32_storage(len); + AutoBuffer buf16_storage(len); + int* buf32 = buf32_storage.data(); + ushort* buf16 = buf16_storage.data(); + + // init from first row + const uchar* src0 = srcmat.ptr(0) + start; + for (int i = 0; i < len; i++) + buf32[i] = (int)src0[i]; + + if (height <= 1) + { + int* dst = dstmat.ptr(0) + start; + for (int i = 0; i < len; i++) + dst[i] = buf32[i]; + return; + } + + memset(buf16, 0, len * sizeof(ushort)); + const int flush_at = 256; // 256*255 = 65280 < 65535 (u16 safe) + int flush_count = 0; + + for (int row = 1; row < height; row++) + { + const uchar* src = srcmat.ptr(row) + start; + int i = 0; + +#if CV_NEON && (defined(__aarch64__) || defined(_M_ARM64)) + // AArch64: vaddw_u8 — single-instruction widening add (u16 += u8) + for (; i <= len - 16; i += 16) + { + uint8x16_t v = vld1q_u8(src + i); + uint16x8_t acc_lo = vld1q_u16(buf16 + i); + uint16x8_t acc_hi = vld1q_u16(buf16 + i + 8); + vst1q_u16(buf16 + i, vaddw_u8(acc_lo, vget_low_u8(v))); + vst1q_u16(buf16 + i + 8, vaddw_high_u8(acc_hi, v)); + } +#else + const int vlanes8 = VTraits::vlanes(); + for (; i <= len - vlanes8; i += vlanes8) + { + v_uint16 lo, hi; + v_expand(vx_load(src + i), lo, hi); + v_store(buf16 + i, v_add(vx_load(buf16 + i), lo)); + v_store(buf16 + i + vlanes16, v_add(vx_load(buf16 + i + vlanes16), hi)); + } +#endif + for (; i < len; i++) + buf16[i] += (ushort)src[i]; + + if (++flush_count >= flush_at || row == height - 1) + { + // flush u16 → u32 + int j = 0; + for (; j <= len - vlanes16; j += vlanes16) + { + v_uint32 lo, hi; + v_expand(vx_load(buf16 + j), lo, hi); + v_store(buf32 + j, v_add(vx_load(buf32 + j), v_reinterpret_as_s32(lo))); + v_store(buf32 + j + vlanes32, v_add(vx_load(buf32 + j + vlanes32), v_reinterpret_as_s32(hi))); + } + for (; j < len; j++) + buf32[j] += (int)buf16[j]; + memset(buf16, 0, len * sizeof(ushort)); + flush_count = 0; + } + } + + int* dst = dstmat.ptr(0) + start; + for (int i = 0; i < len; i++) + dst[i] = buf32[i]; + }; + + parallel_for_(Range(0, width_cn), body, width_cn * CV_ELEM_SIZE(srcmat.depth()) / 64); + v_cleanup(); +} + +// --- uchar → float --- +static void reduceRowSum_8u32f(const Mat& srcmat, Mat& dstmat) +{ + // compute in int, then convert to float + Mat temp(dstmat.rows, dstmat.cols, CV_32SC(srcmat.channels())); + reduceRowSum_8u32s(srcmat, temp); + temp.convertTo(dstmat, dstmat.type()); +} + +// --- 16-bit (ushort/short) → float row reduce --- +template +static void reduceRowSum_16_32f(const Mat& srcmat, Mat& dstmat) +{ + const int width_cn = srcmat.cols * srcmat.channels(); + const int height = srcmat.rows; + const int vlanes16 = VTraits::vlanes(); + const int vlanes32 = vlanes16 / 2; + + auto body = [&](const Range& range) { + const int start = range.start; + const int end = range.end; + const int len = end - start; + float* dst = dstmat.ptr(0) + start; + + int i = 0; + for (; i <= len - vlanes16 * 4; i += vlanes16 * 4) + { + const SrcT* src0 = srcmat.ptr(0) + start + i; + VecDst a0, a1, a2, a3, a4, a5, a6, a7; + v_expand(vx_load(src0), a0, a1); + v_expand(vx_load(src0 + vlanes16), a2, a3); + v_expand(vx_load(src0 + vlanes16*2), a4, a5); + v_expand(vx_load(src0 + vlanes16*3), a6, a7); + for (int row = 1; row < height; row++) + { + const SrcT* src = srcmat.ptr(row) + start + i; + VecDst lo, hi; + v_expand(vx_load(src), lo, hi); + a0 = v_add(a0, lo); a1 = v_add(a1, hi); + v_expand(vx_load(src + vlanes16), lo, hi); + a2 = v_add(a2, lo); a3 = v_add(a3, hi); + v_expand(vx_load(src + vlanes16*2), lo, hi); + a4 = v_add(a4, lo); a5 = v_add(a5, hi); + v_expand(vx_load(src + vlanes16*3), lo, hi); + a6 = v_add(a6, lo); a7 = v_add(a7, hi); + } + v_store(dst + i, v_cvt_f32(v_reinterpret_as_s32(a0))); + v_store(dst + i + vlanes32, v_cvt_f32(v_reinterpret_as_s32(a1))); + v_store(dst + i + vlanes32*2, v_cvt_f32(v_reinterpret_as_s32(a2))); + v_store(dst + i + vlanes32*3, v_cvt_f32(v_reinterpret_as_s32(a3))); + v_store(dst + i + vlanes32*4, v_cvt_f32(v_reinterpret_as_s32(a4))); + v_store(dst + i + vlanes32*5, v_cvt_f32(v_reinterpret_as_s32(a5))); + v_store(dst + i + vlanes32*6, v_cvt_f32(v_reinterpret_as_s32(a6))); + v_store(dst + i + vlanes32*7, v_cvt_f32(v_reinterpret_as_s32(a7))); + } + for (; i <= len - vlanes16; i += vlanes16) + { + VecDst a0, a1; + v_expand(vx_load(srcmat.ptr(0) + start + i), a0, a1); + for (int row = 1; row < height; row++) + { + VecDst lo, hi; + v_expand(vx_load(srcmat.ptr(row) + start + i), lo, hi); + a0 = v_add(a0, lo); a1 = v_add(a1, hi); + } + v_store(dst + i, v_cvt_f32(v_reinterpret_as_s32(a0))); + v_store(dst + i + vlanes32, v_cvt_f32(v_reinterpret_as_s32(a1))); + } + for (; i < len; i++) + { + int val = (int)*(srcmat.ptr(0) + start + i); + for (int row = 1; row < height; row++) + val += (int)*(srcmat.ptr(row) + start + i); + dst[i] = (float)val; + } + }; + + parallel_for_(Range(0, width_cn), body, width_cn * CV_ELEM_SIZE(srcmat.depth()) / 64); + v_cleanup(); +} + +static void reduceRowSum_16u32f(const Mat& srcmat, Mat& dstmat) +{ reduceRowSum_16_32f(srcmat, dstmat); } + +static void reduceRowSum_16s32f(const Mat& srcmat, Mat& dstmat) +{ reduceRowSum_16_32f(srcmat, dstmat); } + +// --- float/double → same type row reduce --- +template +static void reduceRowSum_FP(const Mat& srcmat, Mat& dstmat) +{ + const int width_cn = srcmat.cols * srcmat.channels(); + const int height = srcmat.rows; + const int vlanes = VTraits::vlanes(); + + auto body = [&](const Range& range) { + const int start = range.start; + const int end = range.end; + const int len = end - start; + T* dst = dstmat.ptr(0) + start; + + int i = 0; + for (; i <= len - vlanes * 8; i += vlanes * 8) + { + const T* src0 = srcmat.ptr(0) + start + i; + VecT s0 = vx_load(src0), s1 = vx_load(src0 + vlanes); + VecT s2 = vx_load(src0 + vlanes*2), s3 = vx_load(src0 + vlanes*3); + VecT s4 = vx_load(src0 + vlanes*4), s5 = vx_load(src0 + vlanes*5); + VecT s6 = vx_load(src0 + vlanes*6), s7 = vx_load(src0 + vlanes*7); + for (int row = 1; row < height; row++) + { + const T* src = srcmat.ptr(row) + start + i; + s0 = v_add(s0, vx_load(src)); + s1 = v_add(s1, vx_load(src + vlanes)); + s2 = v_add(s2, vx_load(src + vlanes*2)); + s3 = v_add(s3, vx_load(src + vlanes*3)); + s4 = v_add(s4, vx_load(src + vlanes*4)); + s5 = v_add(s5, vx_load(src + vlanes*5)); + s6 = v_add(s6, vx_load(src + vlanes*6)); + s7 = v_add(s7, vx_load(src + vlanes*7)); + } + v_store(dst + i, s0); + v_store(dst + i + vlanes, s1); + v_store(dst + i + vlanes*2, s2); + v_store(dst + i + vlanes*3, s3); + v_store(dst + i + vlanes*4, s4); + v_store(dst + i + vlanes*5, s5); + v_store(dst + i + vlanes*6, s6); + v_store(dst + i + vlanes*7, s7); + } + for (; i <= len - vlanes; i += vlanes) + { + VecT s0 = vx_load(srcmat.ptr(0) + start + i); + for (int row = 1; row < height; row++) + s0 = v_add(s0, vx_load(srcmat.ptr(row) + start + i)); + v_store(dst + i, s0); + } + for (; i < len; i++) + { + T val = *(srcmat.ptr(0) + start + i); + for (int row = 1; row < height; row++) + val += *(srcmat.ptr(row) + start + i); + dst[i] = val; + } + }; + + parallel_for_(Range(0, width_cn), body, width_cn * CV_ELEM_SIZE(srcmat.depth()) / 64); + v_cleanup(); +} + +static void reduceRowSum_32f32f(const Mat& srcmat, Mat& dstmat) +{ reduceRowSum_FP(srcmat, dstmat); } + +#if CV_SIMD_64F +// --- float → double --- +static void reduceRowSum_32f64f(const Mat& srcmat, Mat& dstmat) +{ + const int width_cn = srcmat.cols * srcmat.channels(); + const int height = srcmat.rows; + const int vlanes32 = VTraits::vlanes(); + const int vlanes64 = VTraits::vlanes(); + + auto body = [&](const Range& range) { + const int start = range.start; + const int end = range.end; + const int len = end - start; + + AutoBuffer buf_storage(len); + double* buf = buf_storage.data(); + + const float* src0 = srcmat.ptr(0) + start; + for (int i = 0; i < len; i++) + buf[i] = (double)src0[i]; + + for (int row = 1; row < height; row++) + { + const float* src = srcmat.ptr(row) + start; + int i = 0; + for (; i <= len - vlanes32; i += vlanes32) + { + v_float32 v = vx_load(src + i); + v_store(buf + i, v_add(vx_load(buf + i), v_cvt_f64(v))); + v_store(buf + i + vlanes64, v_add(vx_load(buf + i + vlanes64), v_cvt_f64_high(v))); + } + for (; i < len; i++) + buf[i] += (double)src[i]; + } + + double* dst = dstmat.ptr(0) + start; + memcpy(dst, buf, len * sizeof(double)); + }; + + parallel_for_(Range(0, width_cn), body, width_cn * CV_ELEM_SIZE(srcmat.depth()) / 64); + v_cleanup(); +} + +static void reduceRowSum_64f64f(const Mat& srcmat, Mat& dstmat) +{ reduceRowSum_FP(srcmat, dstmat); } +#endif // CV_SIMD_64F + +#endif // CV_SIMD || CV_SIMD_SCALABLE + +// ===================================================================== +// Dispatchers +// ===================================================================== + +ReduceSumFunc getReduceCSumFunc(int sdepth, int ddepth) +{ +#if (CV_SIMD || CV_SIMD_SCALABLE) + if (sdepth == CV_8U && ddepth == CV_32S) return reduceColSum_8u32s; + if (sdepth == CV_8U && ddepth == CV_32F) return reduceColSum_8u32f; + if (sdepth == CV_16U && ddepth == CV_32F) return reduceColSum_16u32f; + if (sdepth == CV_16S && ddepth == CV_32F) return reduceColSum_16s32f; + if (sdepth == CV_32F && ddepth == CV_32F) return reduceColSum_32f32f; +#if CV_SIMD_64F + if (sdepth == CV_32F && ddepth == CV_64F) return reduceColSum_32f64f; + if (sdepth == CV_64F && ddepth == CV_64F) return reduceColSum_64f64f; +#endif +#else + CV_UNUSED(sdepth); + CV_UNUSED(ddepth); +#endif + return nullptr; +} + +ReduceSumFunc getReduceRSumFunc(int sdepth, int ddepth) +{ +#if (CV_SIMD || CV_SIMD_SCALABLE) + if (sdepth == CV_8U && ddepth == CV_32S) return reduceRowSum_8u32s; + if (sdepth == CV_8U && ddepth == CV_32F) return reduceRowSum_8u32f; + if (sdepth == CV_16U && ddepth == CV_32F) return reduceRowSum_16u32f; + if (sdepth == CV_16S && ddepth == CV_32F) return reduceRowSum_16s32f; + if (sdepth == CV_32F && ddepth == CV_32F) return reduceRowSum_32f32f; +#if CV_SIMD_64F + if (sdepth == CV_32F && ddepth == CV_64F) return reduceRowSum_32f64f; + if (sdepth == CV_64F && ddepth == CV_64F) return reduceRowSum_64f64f; +#endif +#else + CV_UNUSED(sdepth); + CV_UNUSED(ddepth); +#endif + return nullptr; +} + +#endif // CV_CPU_OPTIMIZATION_DECLARATIONS_ONLY + +CV_CPU_OPTIMIZATION_NAMESPACE_END +} // namespace cv From e458ff16f3c63abd6d90644f700706094bce945d Mon Sep 17 00:00:00 2001 From: Alexander Smorkalov Date: Wed, 22 Apr 2026 09:39:38 +0300 Subject: [PATCH 18/19] Added CI pipeline with vendored 3rdparties and no HAL. --- .github/workflows/PR-4.x.yaml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.github/workflows/PR-4.x.yaml b/.github/workflows/PR-4.x.yaml index 0b25554bfb..7d89a06450 100644 --- a/.github/workflows/PR-4.x.yaml +++ b/.github/workflows/PR-4.x.yaml @@ -12,6 +12,9 @@ jobs: with: workflow_branch: main + Linux-no-HAL: + uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-Linux-NoHAL.yaml@main + Windows: uses: opencv/ci-gha-workflow/.github/workflows/OCV-PR-Windows.yaml@main with: From 273e52ce8d834f6e59f1ed447bf253d475be8ba4 Mon Sep 17 00:00:00 2001 From: Jeevan Mohan Pawar <35501212+jeevan6996@users.noreply.github.com> Date: Wed, 22 Apr 2026 13:42:37 +0100 Subject: [PATCH 19/19] Merge pull request #28818 from jeevan6996:fix-charuco-detectdiamonds-clear-stale-output objdetect: clear stale outputs in CharucoDetector::detectDiamonds #28818 ## Summary - clear `diamondCorners`/`diamondIds` outputs at the beginning of `CharucoDetector::detectDiamonds` - prevent stale detections from previous calls when the current call has fewer than 4 markers or finds no diamonds - add a regression test that pre-fills outputs, calls `detectDiamonds` with 3 markers, and verifies outputs are empty Fixes #28783. ## Testing - built `opencv_test_objdetect` locally - ran `opencv_test_objdetect --gtest_filter=Charuco.detectDiamondsClearsOutputsWithLessThanFourMarkers` - result: PASS --- .../objdetect/src/aruco/charuco_detector.cpp | 30 ++++++++++++++----- .../objdetect/test/test_charucodetection.cpp | 21 +++++++++++++ 2 files changed, 44 insertions(+), 7 deletions(-) diff --git a/modules/objdetect/src/aruco/charuco_detector.cpp b/modules/objdetect/src/aruco/charuco_detector.cpp index 6b19114183..88a567a7a3 100644 --- a/modules/objdetect/src/aruco/charuco_detector.cpp +++ b/modules/objdetect/src/aruco/charuco_detector.cpp @@ -414,7 +414,12 @@ void CharucoDetector::detectDiamonds(InputArray image, OutputArrayOfArrays _diam // stores if the detected markers have been assigned or not to a diamond vector assigned(_markerIds.total(), false); - if(_markerIds.total() < 4ull) return; // a diamond need at least 4 markers + if(_markerIds.total() < 4ull) + { + if (_diamondCorners.needed()) _diamondCorners.release(); + if (_diamondIds.needed()) _diamondIds.release(); + return; // a diamond need at least 4 markers + } // convert input image to grey Mat grey; @@ -518,16 +523,27 @@ void CharucoDetector::detectDiamonds(InputArray image, OutputArrayOfArrays _diam if(diamondIds.size() > 0ull) { // parse output - Mat(diamondIds).copyTo(_diamondIds); + if (_diamondIds.needed()) + { + Mat(diamondIds).copyTo(_diamondIds); + } - _diamondCorners.create((int)diamondCorners.size(), 1, CV_32FC2); - for(unsigned int i = 0; i < diamondCorners.size(); i++) { - _diamondCorners.create(4, 1, CV_32FC2, i, true); - for(int j = 0; j < 4; j++) { - _diamondCorners.getMat(i).at(j) = diamondCorners[i][j]; + if (_diamondCorners.needed()) + { + _diamondCorners.create((int)diamondCorners.size(), 1, CV_32FC2); + for(unsigned int i = 0; i < diamondCorners.size(); i++) { + _diamondCorners.create(4, 1, CV_32FC2, i, true); + for(int j = 0; j < 4; j++) { + _diamondCorners.getMat(i).at(j) = diamondCorners[i][j]; + } } } } + else + { + if (_diamondCorners.needed()) _diamondCorners.release(); + if (_diamondIds.needed()) _diamondIds.release(); + } } void drawDetectedCornersCharuco(InputOutputArray _image, InputArray _charucoCorners, diff --git a/modules/objdetect/test/test_charucodetection.cpp b/modules/objdetect/test/test_charucodetection.cpp index 0e91c69ea9..2d27abcd37 100644 --- a/modules/objdetect/test/test_charucodetection.cpp +++ b/modules/objdetect/test/test_charucodetection.cpp @@ -703,6 +703,27 @@ TEST(Charuco, testmatchImagePoints) } } +TEST(Charuco, detectDiamondsClearsOutputsWithLessThanFourMarkers) +{ + aruco::CharucoBoard board(Size(3, 3), 1.f, 0.5f, aruco::getPredefinedDictionary(aruco::DICT_4X4_50)); + aruco::CharucoDetector detector(board); + + vector> markerCorners = { + {Point2f(10.f, 10.f), Point2f(20.f, 10.f), Point2f(20.f, 20.f), Point2f(10.f, 20.f)}, + {Point2f(30.f, 10.f), Point2f(40.f, 10.f), Point2f(40.f, 20.f), Point2f(30.f, 20.f)}, + {Point2f(10.f, 30.f), Point2f(20.f, 30.f), Point2f(20.f, 40.f), Point2f(10.f, 40.f)} + }; + vector markerIds = {0, 1, 2}; + + vector> diamondCorners = {{Point2f(1.f, 1.f), Point2f(2.f, 1.f), Point2f(2.f, 2.f), Point2f(1.f, 2.f)}}; + vector diamondIds = {Vec4i(0, 1, 2, 3)}; + + detector.detectDiamonds(Mat(), diamondCorners, diamondIds, markerCorners, markerIds); + + EXPECT_TRUE(diamondCorners.empty()); + EXPECT_TRUE(diamondIds.empty()); +} + typedef testing::TestWithParam CharucoDraw; INSTANTIATE_TEST_CASE_P(/**/, CharucoDraw, testing::Values(CV_8UC2, CV_8SC2, CV_16UC2, CV_16SC2, CV_32SC2, CV_32FC2, CV_64FC2)); TEST_P(CharucoDraw, testDrawDetected) {