diff --git a/CMakeLists.txt b/CMakeLists.txt index 22d8c78f26..0edaabe0b0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1465,6 +1465,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) @@ -1493,6 +1495,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/cmake/OpenCVCompilerOptions.cmake b/cmake/OpenCVCompilerOptions.cmake index cb830f7a5d..e43cd9bd79 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() diff --git a/cmake/OpenCVFindLAPACK.cmake b/cmake/OpenCVFindLAPACK.cmake index 4071d501ce..88335b50c9 100644 --- a/cmake/OpenCVFindLAPACK.cmake +++ b/cmake/OpenCVFindLAPACK.cmake @@ -239,11 +239,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" diff --git a/modules/3d/src/precomp.hpp b/modules/3d/src/precomp.hpp index ec4d112654..05d4f7ae99 100755 --- a/modules/3d/src/precomp.hpp +++ b/modules/3d/src/precomp.hpp @@ -47,6 +47,7 @@ # define _SILENCE_CXX17_CODECVT_HEADER_DEPRECATION_WARNING #endif +#include #include "opencv2/core/utility.hpp" #include "opencv2/core/private.hpp" diff --git a/modules/core/CMakeLists.txt b/modules/core/CMakeLists.txt index e7b83358be..68c3237518 100644 --- a/modules/core/CMakeLists.txt +++ b/modules/core/CMakeLists.txt @@ -14,6 +14,7 @@ ocv_add_dispatched_file(minmax SSE2 SSE4_1 AVX2 VSX3 LASX) ocv_add_dispatched_file(nan_mask 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 a703bc5394..9a5ae5fa01 100644 --- a/modules/core/src/matrix_operations.cpp +++ b/modules/core/src/matrix_operations.cpp @@ -327,6 +327,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 { @@ -350,7 +364,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 ) { @@ -801,7 +816,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; @@ -876,7 +894,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/matrix_transform.cpp b/modules/core/src/matrix_transform.cpp index aa1eb29c09..c6eacda3e8 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 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 diff --git a/modules/features/perf/opencl/perf_brute_force_matcher.cpp b/modules/features/perf/opencl/perf_brute_force_matcher.cpp index 439140a5ff..6af06d37ad 100644 --- a/modules/features/perf/opencl/perf_brute_force_matcher.cpp +++ b/modules/features/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/features/src/matchers.cpp b/modules/features/src/matchers.cpp index 249e460e17..2223d59400 100644 --- a/modules/features/src/matchers.cpp +++ b/modules/features/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; @@ -752,6 +756,146 @@ static bool ocl_knnMatch(InputArray query, InputArray _train, std::vector< std:: return false; return true; } + +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::features::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) +{ + 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)) + return false; + + UMat revIdx; + if (!ocl_matchSingle(train, query, revIdx, nullptr, dstType)) + return false; + + 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 +938,15 @@ void BFMatcher::knnMatchImpl( InputArray _queryDescriptors, std::vector 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 diff --git a/modules/features/test/test_matchers_algorithmic.cpp b/modules/features/test/test_matchers_algorithmic.cpp index 0679d9ce38..317b15beee 100644 --- a/modules/features/test/test_matchers_algorithmic.cpp +++ b/modules/features/test/test_matchers_algorithmic.cpp @@ -629,4 +629,47 @@ TEST(Features2d_DMatch, issue_17771) EXPECT_NO_THROW(ubf->knnMatch(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 diff --git a/modules/imgcodecs/CMakeLists.txt b/modules/imgcodecs/CMakeLists.txt index 68841e93b6..56ddc96f6d 100644 --- a/modules/imgcodecs/CMakeLists.txt +++ b/modules/imgcodecs/CMakeLists.txt @@ -194,4 +194,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/include/opencv2/imgcodecs.hpp b/modules/imgcodecs/include/opencv2/imgcodecs.hpp index ac5591c8c7..21b4c7461c 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. @@ -541,8 +541,11 @@ can be saved using this function, with these exceptions: 64-bit unsigned (CV_64U), 64-bit signed (CV_64S), 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_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/src/grfmt_tiff.cpp b/modules/imgcodecs/src/grfmt_tiff.cpp index 06a3f3c3b0..31f48d72a0 100644 --- a/modules/imgcodecs/src/grfmt_tiff.cpp +++ b/modules/imgcodecs/src/grfmt_tiff.cpp @@ -1323,8 +1323,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)) @@ -1431,14 +1435,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; @@ -1493,7 +1520,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; } @@ -1513,10 +1539,10 @@ bool TiffEncoder::writeLibTiff( const std::vector& img_vec, const std::vect // Predictor 2 for 64-bit is supported at v4.4.0 or later. // See https://libtiff.gitlab.io/libtiff/releases/v4.4.0.html #if TIFFLIB_VERSION < 20220520 /* Magic number of libtiff v4.4.0 */ - if ( (bitsPerChannel == 64) && (predictor == PREDICTOR_HORIZONTAL /* 2 */) ) + if ( (bitsPerChannel == 64) && (page_predictor == PREDICTOR_HORIZONTAL /* 2 */) ) { CV_LOG_ONCE_WARNING(NULL, "Predictor 2(HORIZONTAL) for 64-bit is supported at v4.4.0 or later, so it is fallbacked to 0(NONE)"); - predictor = PREDICTOR_NONE; + page_predictor = PREDICTOR_NONE; } #endif @@ -1541,7 +1567,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) @@ -1626,6 +1652,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_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); diff --git a/modules/imgcodecs/test/test_tiff.cpp b/modules/imgcodecs/test/test_tiff.cpp index 4850a22dfa..8afd5ec5f2 100644 --- a/modules/imgcodecs/test/test_tiff.cpp +++ b/modules/imgcodecs/test/test_tiff.cpp @@ -937,7 +937,7 @@ Imgcodes_Tiff_TypeAndComp all_types[] = { { CV_32SC1, true }, { CV_32SC3, true }, { CV_32SC4, true }, { CV_64UC1, true }, { CV_64UC3, true }, { CV_64UC4, true }, { CV_64SC1, true }, { CV_64SC3, true }, { CV_64SC4, 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 }; @@ -1296,6 +1296,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 diff --git a/modules/imgproc/src/shapedescr.cpp b/modules/imgproc/src/shapedescr.cpp index 302a718862..764b0a6d7a 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); 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; diff --git a/modules/objdetect/src/aruco/charuco_detector.cpp b/modules/objdetect/src/aruco/charuco_detector.cpp index d2accdf011..4f733403db 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; @@ -522,16 +527,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 793d4dcfa1..2dfbcb199f 100644 --- a/modules/objdetect/test/test_charucodetection.cpp +++ b/modules/objdetect/test/test_charucodetection.cpp @@ -705,6 +705,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) { 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); } diff --git a/modules/video/include/opencv2/video/tracking.hpp b/modules/video/include/opencv2/video/tracking.hpp index 542f02a085..acd3eb753c 100644 --- a/modules/video/include/opencv2/video/tracking.hpp +++ b/modules/video/include/opencv2/video/tracking.hpp @@ -428,6 +428,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/snippets/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 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/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 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