Merge branch 4.x

2026-07-21 19:33:03 +04:00 · 2026-04-23 11:00:34 +03:00
parent 323ea1782d 273e52ce8d
commit 7bde3f9ae9
38 changed files with 3437 additions and 293 deletions
@@ -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)
@@ -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()
@@ -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"
@@ -47,6 +47,7 @@
 #  define _SILENCE_CXX17_CODECVT_HEADER_DEPRECATION_WARNING
 #endif

+#include <opencv2/core/utils/logger.hpp>
 #include "opencv2/core/utility.hpp"

 #include "opencv2/core/private.hpp"
@@ -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
@@ -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 <typename T, typename WT, typename Op>
+struct ReduceR_SIMD
+{
+    int operator()(const T*, int start, int, WT*, const Op&) const
+    {
+        return start;
+    }
+};
+
+
 template<typename T, typename ST, typename WT, class Op, class OpInit>
 class ReduceR_Invoker : public ParallelLoopBody
 {
@@ -350,7 +364,8 @@ public:
    for( ; --height; )
    {
        src += srcstep;
-        i = range.start;
+        ReduceR_SIMD<T, WT, Op> 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;
@@ -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
@@ -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
@@ -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
@@ -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<DMatch> 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

@@ -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<DescriptorMatcher> BFMatcher::clone( bool emptyTrainData ) const
 static bool ocl_match(InputArray query, InputArray _train, std::vector< std::vector<DMatch> > &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<DMatch> >& 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<int>();
+    const float* dist = fwdDistCPU.ptr<float>();
+    const int* revBestData = revBestCPU.ptr<int>();
+
+    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<DMatch>(1, DMatch(q, t, 0, dist[q])));
+            }
+        }
+    }
+    return true;
+}
+
+static bool ocl_matchWithCrossCheck(InputArray query, InputArray train,
+        std::vector< std::vector<DMatch> >& 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<int>();
+    const int* rev  = revIdxCPU.ptr<int>();
+    const float* dist = fwdDistCPU.ptr<float>();
+
+    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<DMatch>(1, DMatch(q, t, 0, dist[q])));
+        }
+    }
+    return true;
+}
 #endif

 void BFMatcher::knnMatchImpl( InputArray _queryDescriptors, std::vector<std::vector<DMatch> >& matches, int knn,
@@ -794,21 +938,15 @@ void BFMatcher::knnMatchImpl( InputArray _queryDescriptors, std::vector<std::vec
    {
        if(knn == 1)
        {
-            if(trainDescCollection.empty())
+            InputArray train = trainDescCollection.empty() ?
+                (InputArray)utrainDescCollection[0] : (InputArray)trainDescCollection[0];
+            bool oclOk = crossCheck ?
+                ocl_matchWithCrossCheck(_queryDescriptors, train, matches, normType) :
+                ocl_match(_queryDescriptors, train, matches, normType);
+            if (oclOk)
            {
-                if(ocl_match(_queryDescriptors, utrainDescCollection[0], matches, normType))
-                {
-                    CV_IMPL_ADD(CV_IMPL_OCL);
-                    return;
-                }
-            }
-            else
-            {
-                if(ocl_match(_queryDescriptors, trainDescCollection[0], matches, normType))
-                {
-                    CV_IMPL_ADD(CV_IMPL_OCL);
-                    return;
-                }
+                CV_IMPL_ADD(CV_IMPL_OCL);
+                return;
            }
        }
        else
@@ -557,4 +557,138 @@ __kernel void BruteForceMatch_knnMatch(
        bestTrainIdx[queryIdx] = (int2)(myBestTrainIdx1, myBestTrainIdx2);
        bestDistance[queryIdx] = (float2)(myBestDistance1, myBestDistance2);
    }
-}
+}
+
+#ifdef HAVE_INT64_ATOMICS
+#pragma OPENCL EXTENSION cl_khr_int64_base_atomics : enable
+
+__kernel void BruteForceMatch_CrossCheckMatch(
+    __global T *query,
+    __global T *train,
+    __global int *bestTrainIdx,
+    __global float *bestDistance,
+    __global ulong *revBest,
+    int query_rows,
+    int query_cols,
+    int train_rows,
+    int train_cols,
+    int step
+)
+{
+    const int lidx = get_local_id(0);
+    const int lidy = get_local_id(1);
+    const int groupidx = get_group_id(0);
+
+    const int queryIdx = mad24(BLOCK_SIZE, groupidx, lidy);
+    const int queryOffset = min(queryIdx, query_rows - 1) * step;
+    __global TN *query_vec = (__global TN *)(query + queryOffset);
+    query_cols /= kercn;
+
+    __local float sharebuffer[SHARED_MEM_SZ];
+    __local value_type *s_query = (__local value_type *)sharebuffer;
+
+#if 0 < MAX_DESC_LEN
+    __local value_type *s_train = (__local value_type *)sharebuffer + BLOCK_SIZE * MAX_DESC_LEN;
+    #pragma unroll
+    for (int i = 0; i < MAX_DESC_LEN / BLOCK_SIZE; i++)
+    {
+        const int loadx = mad24(BLOCK_SIZE, i, lidx);
+        s_query[mad24(MAX_DESC_LEN, lidy, loadx)] = loadx < query_cols ? query_vec[loadx] : 0;
+    }
+#else
+    __local value_type *s_train = (__local value_type *)sharebuffer + BLOCK_SIZE_ODD * BLOCK_SIZE;
+    const int s_query_i = mad24(BLOCK_SIZE_ODD, lidy, lidx);
+    const int s_train_i = mad24(BLOCK_SIZE_ODD, lidx, lidy);
+#endif
+
+    float myBestDistance = MAX_FLOAT;
+    int myBestTrainIdx = -1;
+
+    for (int t = 0, endt = (train_rows + BLOCK_SIZE - 1) / BLOCK_SIZE; t < endt; t++)
+    {
+        result_type result = 0;
+
+        const int trainOffset = min(mad24(BLOCK_SIZE, t, lidy), train_rows - 1) * step;
+        __global TN *train_vec = (__global TN *)(train + trainOffset);
+#if 0 < MAX_DESC_LEN
+        #pragma unroll
+        for (int i = 0; i < MAX_DESC_LEN / BLOCK_SIZE; i++)
+        {
+            const int loadx = mad24(BLOCK_SIZE, i, lidx);
+            s_train[mad24(BLOCK_SIZE, lidx, lidy)] = loadx < train_cols ? train_vec[loadx] : 0;
+
+            barrier(CLK_LOCAL_MEM_FENCE);
+
+            result += reduce_multi_block(s_query, s_train, i, lidx, lidy);
+
+            barrier(CLK_LOCAL_MEM_FENCE);
+        }
+#else
+        for (int i = 0, endq = (query_cols + BLOCK_SIZE - 1) / BLOCK_SIZE; i < endq; i++)
+        {
+            const int loadx = mad24(BLOCK_SIZE, i, lidx);
+            if (loadx < query_cols)
+            {
+                s_query[s_query_i] = query_vec[loadx];
+                s_train[s_train_i] = train_vec[loadx];
+            }
+            else
+            {
+                s_query[s_query_i] = 0;
+                s_train[s_train_i] = 0;
+            }
+
+            barrier(CLK_LOCAL_MEM_FENCE);
+
+            result += reduce_block_match(s_query, s_train, lidx, lidy);
+
+            barrier(CLK_LOCAL_MEM_FENCE);
+        }
+#endif
+        result = DIST_RES(result);
+
+        const int trainIdx = mad24(BLOCK_SIZE, t, lidx);
+
+        if (queryIdx < query_rows && trainIdx < train_rows)
+        {
+            if (result < myBestDistance)
+            {
+                myBestDistance = result;
+                myBestTrainIdx = trainIdx;
+            }
+
+            uint dist_bits = as_uint(result);
+            ulong packed = ((ulong)dist_bits << 32) | (ulong)queryIdx;
+            atom_min(&revBest[trainIdx], packed);
+        }
+    }
+
+    barrier(CLK_LOCAL_MEM_FENCE);
+
+    __local float *s_distance = (__local float *)sharebuffer;
+    __local int *s_trainIdx = (__local int *)(sharebuffer + BLOCK_SIZE_ODD * BLOCK_SIZE);
+
+    s_distance += lidy * BLOCK_SIZE_ODD;
+    s_trainIdx += lidy * BLOCK_SIZE_ODD;
+    s_distance[lidx] = myBestDistance;
+    s_trainIdx[lidx] = myBestTrainIdx;
+
+    barrier(CLK_LOCAL_MEM_FENCE);
+
+    #pragma unroll
+    for (int k = 0; k < BLOCK_SIZE; k++)
+    {
+        if (myBestDistance > s_distance[k])
+        {
+            myBestDistance = s_distance[k];
+            myBestTrainIdx = s_trainIdx[k];
+        }
+    }
+
+    if (queryIdx < query_rows && lidx == 0)
+    {
+        bestTrainIdx[queryIdx] = myBestTrainIdx;
+        bestDistance[queryIdx] = myBestDistance;
+    }
+}
+#endif
@@ -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<BFMatcher> cpuMatcher = BFMatcher::create(NORM_L2, true /*crossCheck*/);
+    vector<DMatch> 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<BFMatcher> oclMatcher = BFMatcher::create(NORM_L2, true /*crossCheck*/);
+    vector<DMatch> 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
@@ -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()
@@ -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.
@@ -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<uint8_t>(prim_id),
+                    static_cast<uint8_t>(tran_id),
+                    static_cast<uint8_t>(matrix_id),
+                    static_cast<uint8_t>(video_full_range_flag)
+                };
+                auto& out = m_metadata[IMAGE_METADATA_CICP];
+                out.insert(out.end(), cicp_data, cicp_data + 4);
+            }
+#endif

            result = true;
        }
@@ -1323,8 +1323,12 @@ bool TiffEncoder::writeLibTiff( const std::vector<Mat>& img_vec, const std::vect
    cv::Ptr<void> 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<Mat>& 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<Mat>& 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<Mat>& 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<Mat>& 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<float>(i), strip_size * sizeof(float)) != (tsize_t)-1);
+    }
+    CV_TIFF_CHECK_CALL(TIFFWriteDirectory(tif));
+    return true;
+}
+
 bool TiffEncoder::writemulti(const std::vector<Mat>& img_vec, const std::vector<int>& params)
 {
    return writeLibTiff(img_vec, params);
@@ -134,6 +134,7 @@ public:
 protected:
    bool writeLibTiff( const std::vector<Mat>& img_vec, const std::vector<int>& 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
@@ -144,6 +144,17 @@ namespace opencv_test { namespace {
        return iccp_data;
    }

+#ifdef OPENCV_IMGCODECS_PNG_WITH_cICP
+    static std::vector<uchar> 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<int> 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());
 }

@@ -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);
@@ -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_Values> 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<int> params;
+    if (compression > 0)
+    {
+      params.push_back(IMWRITE_TIFF_COMPRESSION);
+      params.push_back(compression);
+    }
+
+    std::vector<unsigned char> 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
@@ -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<double, 3, 1> 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<double>(0,0) * eVec.at<double>(0,2) - eVec.at<double>(0,1) * eVec.at<double>(0,1));
        cond[1]=(4.0 * eVec.at<double>(1,0) * eVec.at<double>(1,2) - eVec.at<double>(1,1) * eVec.at<double>(1,1));
        cond[2]=(4.0 * eVec.at<double>(2,0) * eVec.at<double>(2,2) - eVec.at<double>(2,1) * eVec.at<double>(2,1));
@@ -793,6 +793,18 @@ cv::RotatedRect cv::fitEllipseDirect( InputArray _points )
        } else {
            i = (cond[0]<cond[2]) ? 2 : 0;
        }
+        // Check that the ellipse condition 4ac-b^2 is meaningfully positive
+        // (not just positive due to floating-point noise from complex eigenvalues)
+        {
+            double v0 = eVec.at<double>(i,0), v1 = eVec.at<double>(i,1), v2 = eVec.at<double>(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<double>(i,0)*eVec.at<double>(i,0) + eVec.at<double>(i,1)*eVec.at<double>(i,1) + eVec.at<double>(i,2)*eVec.at<double>(i,2));
        if (((eVec.at<double>(i,0)<0.0  ? -1 : 1) * (eVec.at<double>(i,1)<0.0  ? -1 : 1) * (eVec.at<double>(i,2)<0.0  ? -1 : 1)) <= 0.0) {
                norm=-1.0*norm;
@@ -337,6 +337,26 @@ TEST(Imgproc_FitEllipseAMS_Issue_7, accuracy) {
    EXPECT_TRUE(checkEllipse(ellipseAMSTest, ellipseAMSTrue, tol));
 }

+TEST(Imgproc_FitEllipseAMS_NearCircular, accuracy)
+{
+    std::vector<cv::Point2f> 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<Point2f> pts({{-300, 100}, {-200, 100}, {-100, 100}, {0, 100}, {100, 100}, {200, 100}, {300, 100}});
    const RotatedRect el = fitEllipseAMS(pts);
@@ -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<cv::Point2f> 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<Point2f> pts({{-300, 100}, {-200, 100}, {-100, 100}, {0, 100}, {100, 100}, {200, 100}, {300, 100}});
    const RotatedRect el = fitEllipseDirect(pts);
@@ -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<unsigned char>(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;
@@ -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<bool> 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<Point2f>(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<Point2f>(j) = diamondCorners[i][j];
+                }
            }
        }
    }
+    else
+    {
+        if (_diamondCorners.needed()) _diamondCorners.release();
+        if (_diamondIds.needed()) _diamondIds.release();
+    }
 }

 void drawDetectedCornersCharuco(InputOutputArray _image, InputArray _charucoCorners,
@@ -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<vector<Point2f>> 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<int> markerIds = {0, 1, 2};
+
+    vector<vector<Point2f>> diamondCorners = {{Point2f(1.f, 1.f), Point2f(2.f, 1.f), Point2f(2.f, 2.f), Point2f(1.f, 2.f)}};
+    vector<Vec4i> diamondIds = {Vec4i(0, 1, 2, 3)};
+
+    detector.detectDiamonds(Mat(), diamondCorners, diamondIds, markerCorners, markerIds);
+
+    EXPECT_TRUE(diamondCorners.empty());
+    EXPECT_TRUE(diamondIds.empty());
+}
+
 typedef testing::TestWithParam<int> 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) {
@@ -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);
 }

@@ -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<int> itersPerLevel = std::vector<int>();
+    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
 */
@@ -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<MotionType, ReadFlag> TestParams;
+typedef std::tuple<MotionType, MultiScaleFlag> TestParamsMS;
 typedef perf::TestBaseWithParam<TestParams> ECCPerfTest;
+typedef perf::TestBaseWithParam<TestParamsMS> ECCPerfTestMS;
+
+typedef std::tuple<MotionType, ReadFlag> 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_<float>(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_<float>(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_<float>(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_<float>(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
@@ -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<cv::Mat> MatPyramid;
+
+template<int motionType> struct MotionTraits {};
+
+template<> struct MotionTraits<MOTION_TRANSLATION> {
+    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<typename elemtype>
+    static constexpr std::array<float, paramAmount> 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<float, paramAmount>{GX, GY};
+#undef GX
+#undef GY
+    }
+};
+
+template<> struct MotionTraits<MOTION_EUCLIDEAN> {
+    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<typename elemtype>
+    static constexpr std::array<float, paramAmount> 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<float, paramAmount>{GZ, GX, GY};
+#undef GX
+#undef GY
+#undef HATX
+#undef HATY
+#undef GZ
+    }
+};
+
+template<> struct MotionTraits<MOTION_AFFINE> {
+    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<typename elemtype>
+    static constexpr std::array<float, paramAmount> 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<float, paramAmount>{GX * col, GY * col, GX * row, GY * row, GX, GY};
+#undef GX
+#undef GY
+    }
+};
+
+template<> struct MotionTraits<MOTION_HOMOGRAPHY> {
+    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<typename elemtype>
+    static constexpr std::array<float, paramAmount> 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<float, paramAmount>{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<int N, class F>
+class constexprForClass
+{
+public:
+    static inline void execute(F&& fVal) {
+        constexprForClass<N-1, F>::execute(std::forward<F>(fVal));
+        fVal(N-1);
+    }
+};
+
+template<class F>
+class constexprForClass<0, F>
+{
+public:
+    static inline void execute(F&&) {}
+};
+
+template<int N, class F>
+void constexprFor(F&& fVal) {
+    constexprForClass<N, F>::execute(std::forward<F>(fVal));
+}
+template<int R, int C, class F>
+class constexprForUpperTriangleClassOneRow
+{
+public:
+    static inline void execute(F&& fVal) {
+        constexprForUpperTriangleClassOneRow<R, C-1, F>::execute(std::forward<F>(fVal));
+        fVal(R, R + C - 1);
+    }
+};
+
+template<int R, class F>
+class constexprForUpperTriangleClassOneRow<R, 0, F>
+{
+public:
+    static inline void execute(F&&) {}
+};
+
+template<int R, int D, class F>
+class constexprForUpperTriangleClass
+{
+public:
+    static inline void execute(F&& fVal) {
+        constexprForUpperTriangleClass<R-1, D, F>::execute(std::forward<F>(fVal));
+        constexprForUpperTriangleClassOneRow<R-1, D-R+1, F>::execute(std::forward<F>(fVal));
+    }
+};
+
+template<int D, class F>
+class constexprForUpperTriangleClass<0, D, F>
+{
+public:
+    static inline void execute(F&&) {}
+};
+
+template<int M, class F>
+void constexprForUpperTriangle(F&& fVal) {
+    constexprForUpperTriangleClass<M,M,F>::execute(std::forward<F>(fVal));
+}
+
+template<int MotionType>
+constexpr int hessianRowStart(int row) {
+    return row == 0 ? 0 : (MotionTraits<MotionType>::paramAmount - row + 1 + hessianRowStart<MotionType>(row - 1));
+}
+
+template<int motionType, typename elemtype>
+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<float, elemtype>::value, "imageHessianProjECC: f16 is not supported yet");
+#define HESSIAN_PARAMS (MotionTraits<motionType>::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<float, elemtype>::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<Matx<double, MotionTraits<motionType>::paramAmount, MotionTraits<motionType>::paramAmount> > hessPs(stripesAmount);
+    std::vector<Vec<double, MotionTraits<motionType>::paramAmount> > iprojs(stripesAmount);
+    std::vector<Vec<double, MotionTraits<motionType>::paramAmount> > tprojs(stripesAmount);
+    std::vector<Vec<double, MotionTraits<motionType>::paramAmount> > projSubs(stripesAmount);
+    std::vector<double> correlations(stripesAmount, 0.);
+
+    std::vector<double> sampSums(stripesAmount, 0);
+    std::vector<double> sampSqSums(stripesAmount, 0);
+    std::vector<double> refSums(stripesAmount, 0);
+    std::vector<double> refSqSums(stripesAmount, 0);
+    std::vector<int> nzs(stripesAmount, 0);
+    std::vector<double> sampMaskedSums(stripesAmount, 0);
+    std::vector<double> refMaskedSums(stripesAmount, 0);
+
+    double a00 = map.at<double>(0, 0);
+    double a01 = map.at<double>(0, 1);
+    double a02 = map.at<double>(0, 2);
+    double a10 = map.at<double>(1, 0);
+    double a11 = map.at<double>(1, 1);
+    double a12 = map.at<double>(1, 2);
+    double a20 = 0;
+    double a21 = 0;
+    double a22 = 0;
+    if (motionType == MOTION_HOMOGRAPHY) {
+        a20 = map.at<double>(2, 0);
+        a21 = map.at<double>(2, 1);
+        a22 = map.at<double>(2, 2);
+    }
+
+    const elemtype* samplePtr0 = sampleWithGrad.ptr<elemtype>(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<elemtype>(y);
+
+            std::array<float, (HESSIAN_PARAMS * HESSIAN_PARAMS + HESSIAN_PARAMS) / 2> hessPcache{};
+            std::array<float, HESSIAN_PARAMS> iprojCache{};
+            std::array<float, HESSIAN_PARAMS> tprojCache{};
+            std::array<float, HESSIAN_PARAMS> 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<motionType>::tailHandlerGetCoord(sx, sy, denominator, x, numeratorX0, numeratorY0,
+                                                denominator0, (float)a00, (float)a10, (float)a20);
+                const unsigned int x0 = (interpolation == INTER_LINEAR) ? static_cast<int>(std::floor(sx)) : saturate_cast<unsigned>(sx);
+                const unsigned int y0 = (interpolation == INTER_LINEAR) ? static_cast<int>(std::floor(sy)) : saturate_cast<unsigned>(sy);
+                if(interpolation == INTER_LINEAR && (static_cast<int>(x0 < xcond) & static_cast<int>(y0 < ycond)) == 0)
+                    continue;
+                if (interpolation == INTER_NEAREST && (static_cast<int>(x0 < xcond) & static_cast<int>(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<float, HESSIAN_PARAMS> jac = MotionTraits<motionType>::fillJacobian(x, y, sx, sy,
+                                                                                        fVal, gx,
+                                                                                        gy, (float)a00,
+                                                                                        (float)a10, denominator);
+                constexprForUpperTriangle<HESSIAN_PARAMS>([&](int row_i, int col_i) {
+                    hessPcache[hessianRowStart<motionType>(row_i) + (col_i - row_i)] += jac[row_i] * jac[col_i];
+                });
+                constexprFor<HESSIAN_PARAMS>([&](int elem) {
+                    iprojCache[elem] += jac[elem] * sampleVal;
+                    tprojCache[elem] += jac[elem] * refVal;
+                    projSubCache[elem] += jac[elem] * fVal;
+                });
+                correlations[stripeIdx] += sampleVal * refVal;
+            }
+
+            constexprForUpperTriangle<HESSIAN_PARAMS>([&](int row, int col) {
+                hessPs[stripeIdx](row, col) += hessPcache[hessianRowStart<motionType>(row) + (col - row)];
+            });
+            constexprFor<HESSIAN_PARAMS>([&](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<double>(0);
+    double* sampleProjPtr = sampleProj.ptr<double>(0);
+    double* refProjPtr = refProj.ptr<double>(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<HESSIAN_PARAMS>([&](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<double>(0);
+    const double* updatePtr = update.ptr<double>(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<MOTION_TRANSLATION, float>(map,
+                                                                       sampleWithGrad,
+                                                                       reference,
+                                                                       sampSum,
+                                                                       sampSqSum,
+                                                                       referenceSum,
+                                                                       referenceSqSum,
+                                                                       nz,
+                                                                       hessian,
+                                                                       sampleProjection,
+                                                                       referenceProjection,
+                                                                       deltaY,
+                                                                       interpolation);
+        } else if (motionType == MOTION_EUCLIDEAN) {
+            correlation = imageHessianProjECC<MOTION_EUCLIDEAN, float>(map,
+                                                                       sampleWithGrad,
+                                                                       reference,
+                                                                       sampSum,
+                                                                       sampSqSum,
+                                                                       referenceSum,
+                                                                       referenceSqSum,
+                                                                       nz,
+                                                                       hessian,
+                                                                       sampleProjection,
+                                                                       referenceProjection,
+                                                                       deltaY,
+                                                                       interpolation);
+        } else if (motionType == MOTION_AFFINE) {
+            correlation = imageHessianProjECC<MOTION_AFFINE, float>(map,
+                                                                       sampleWithGrad,
+                                                                       reference,
+                                                                       sampSum,
+                                                                       sampSqSum,
+                                                                       referenceSum,
+                                                                       referenceSqSum,
+                                                                       nz,
+                                                                       hessian,
+                                                                       sampleProjection,
+                                                                       referenceProjection,
+                                                                       deltaY,
+                                                                       interpolation);
+        } else {
+            correlation = imageHessianProjECC<MOTION_HOMOGRAPHY, float>(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<float>();
+        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<float>(row);
+                const float* samplePrevLine = sample.ptr<float>(std::max(row - 1, 0));
+                const float* sampleNextLine = sample.ptr<float>(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<float>(row);
+                const float* srcPtr = image.ptr<float>(row);
+                const uint8_t* maskPtr = !mask.empty() ? mask.ptr<uint8_t>(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<double>(0, 0) = 1.f / scale;
+        invertScaleMat.at<double>(1, 1) = 1.f / scale;
+        invertScaleMat.at<double>(2, 2) = 1.f;
+        Mat scaleMatrix = invertScaleMat.clone();
+        scaleMatrix.at<double>(0, 0) = scale;
+        scaleMatrix.at<double>(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<double*>(warpMatrix.data))[mel] /=
+                (reinterpret_cast<double*>(warpMatrix.data))[8];
+        }
+        (reinterpret_cast<double*>(warpMatrix.data))[8] = 1.f;
+    } else {
+        warpMatrix.at<double>(0, 2) *= scale;
+        warpMatrix.at<double>(1, 2) *= scale;
+    }
+}
+
+static void checkParams(const MatPyramid& referencePyramid,
+                 const MatPyramid& samplePyramid,
+                 Mat& map,
+                 std::vector<int>& 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<int>(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<int>(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<int> 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<MOTION_TRANSLATION>::paramAmount;
+            break;
+        case MOTION_EUCLIDEAN:
+            nparams = MotionTraits<MOTION_EUCLIDEAN>::paramAmount;
+            break;
+        case MOTION_AFFINE:
+            nparams = MotionTraits<MOTION_AFFINE>::paramAmount;
+            break;
+        case MOTION_HOMOGRAPHY:
+            nparams = MotionTraits<MOTION_HOMOGRAPHY>::paramAmount;
+            break;
+        default:
+            CV_Error(Error::StsBadArg, "Incorrect motion type");
+    }
+
+    const std::vector<int> numberOfIterations = ((eccParams.criteria.type & TermCriteria::COUNT) != 0)
+                                                    ? itersPerLevelCopy
+                                                    : std::vector<int>(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. */
@@ -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_<float>(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_<float>(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_<float>(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_<float>(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_<float>(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_<float>(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_<float>(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_<float>(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_<float>(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_<float>(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_<float>(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_<float>(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_<float>(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_<float>(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
@@ -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
+#pragma clang diagnostic pop
@@ -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);
@@ -1105,7 +1105,7 @@ bool CvCapture_FFMPEG::open(const char* _filename, int index, const Ptr<IStreamR
                {
                    enableAlpha = false;
                }
-                if (value == CV_8UC4)
+                else if (value == CV_8UC4)
                {
                    enableAlpha = true;
                }
@@ -845,6 +845,25 @@ TEST(videoio_ffmpeg, create_with_property_badarg)
    EXPECT_FALSE(cap.isOpened());
 }

+// requires FFmpeg wrapper rebuild on Windows
+#ifndef _WIN32
+TEST(videoio_ffmpeg, open_with_format_cv8uc3)
+{
+    if (!videoio_registry::hasBackend(CAP_FFMPEG))
+        throw SkipTestException("FFmpeg backend was not found");
+
+    string video_file = findDataFile("video/big_buck_bunny.mp4");
+    VideoCapture cap(video_file, CAP_FFMPEG, {
+        CAP_PROP_FORMAT, CV_8UC3
+    });
+    ASSERT_TRUE(cap.isOpened());
+    EXPECT_EQ(cap.get(CAP_PROP_FORMAT), CV_8UC3);
+    Mat frame;
+    ASSERT_TRUE(cap.read(frame));
+    EXPECT_EQ(frame.channels(), 3);
+}
+#endif
+
 // related issue: https://github.com/opencv/opencv/issues/16821
 TEST(videoio_ffmpeg, DISABLED_open_from_web)
 {
@@ -1236,4 +1236,43 @@ inline static std::string stream_capture_ffmpeg_name_printer(const testing::Test

 INSTANTIATE_TEST_CASE_P(videoio, stream_capture_ffmpeg, testing::Values("h264", "h265", "mjpg.avi"), stream_capture_ffmpeg_name_printer);

+// Seek by frame index must be accurate on non-integer frame rates.
+typedef testing::TestWithParam<VideoCaptureAPIs> 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