From 9f101a126f1c7b0f71f56491ce751e009fce8fff Mon Sep 17 00:00:00 2001 From: 4ekmah Date: Mon, 20 Apr 2026 17:22:19 +0300 Subject: [PATCH] Merge pull request #28802 from 4ekmah:pyr_ecc Multiscale ECC #28802 OpenCV Extra: https://github.com/opencv/opencv_extra/pull/1338 ### Pull Request Readiness Checklist See details at https://github.com/opencv/opencv/wiki/How_to_contribute#making-a-good-pull-request - [x] I agree to contribute to the project under Apache 2 License. - [x] To the best of my knowledge, the proposed patch is not based on a code under GPL or another license that is incompatible with OpenCV - [x] The PR is proposed to the proper branch - [ ] There is a reference to the original bug report and related work - [x] There is accuracy test, performance test and test data in opencv_extra repository, if applicable Patch to opencv_extra has the same branch name. - [x] The feature is well documented and sample code can be built with the project CMake --- .../video/include/opencv2/video/tracking.hpp | 82 ++ modules/video/perf/perf_ecc.cpp | 60 ++ modules/video/src/eccms.cpp | 885 ++++++++++++++++++ modules/video/test/test_ecc.cpp | 368 ++++---- 4 files changed, 1211 insertions(+), 184 deletions(-) create mode 100644 modules/video/src/eccms.cpp diff --git a/modules/video/include/opencv2/video/tracking.hpp b/modules/video/include/opencv2/video/tracking.hpp index 2ed27f505d..f62660494c 100644 --- a/modules/video/include/opencv2/video/tracking.hpp +++ b/modules/video/include/opencv2/video/tracking.hpp @@ -419,6 +419,88 @@ CV_EXPORTS_W double findTransformECCWithMask( InputArray templateImage, TermCriteria criteria = TermCriteria(TermCriteria::COUNT + TermCriteria::EPS, 50, 1e-6), int gaussFiltSize = 5 ); +/** @brief struct ECCParameters is used by findTransformECCMultiScale + +@param motionType parameter, specifying the type of motion: + - **MOTION_TRANSLATION** sets a translational motion model; warpMatrix is \f$2\times 3\f$ with + the first \f$2\times 2\f$ part being the unity matrix and the rest two parameters being + estimated. + - **MOTION_EUCLIDEAN** sets a Euclidean (rigid) transformation as motion model; three + parameters are estimated; warpMatrix is \f$2\times 3\f$. + - **MOTION_AFFINE** sets an affine motion model (DEFAULT); six parameters are estimated; + warpMatrix is \f$2\times 3\f$. + - **MOTION_HOMOGRAPHY** sets a homography as a motion model; eight parameters are + estimated;\`warpMatrix\` is \f$3\times 3\f$. +@param criteria parameter, specifying the termination criteria of the ECC algorithm; +criteria.epsilon defines the threshold of the increment in the correlation coefficient between two +iterations (a negative criteria.epsilon makes criteria.maxcount the only termination criterion). +Default values are shown in the declaration above. +@param itersPerLevel Criterion extension: distribution of iterations limit over pyramid levels. +Can be empty, in this case, this algorithm will use criteria.maxCount on each level. +@param gaussFiltSize An optional value indicating size of gaussian blur filter; (DEFAULT: 5) +@param nlevels An optional value indicating amount of levels in the pyramid; (DEFAULT: 4) +@param interpolation Type of warp interpolation. Possible values are INTER_NEAREST and INTER_LINEAR. +Affects accuracy, especially when motionType == MOTION_TRANSLATION. (DEFAULT: INTER_LINEAR) + */ +struct CV_EXPORTS_W_SIMPLE ECCParameters +{ + CV_WRAP ECCParameters() {} + CV_PROP_RW int motionType = MOTION_AFFINE; + CV_PROP_RW cv::TermCriteria criteria = TermCriteria(TermCriteria::COUNT + TermCriteria::EPS, 50, 1e-6); + CV_PROP_RW std::vector itersPerLevel = std::vector(); + CV_PROP_RW int gaussFiltSize = 5; + CV_PROP_RW int nlevels = 4; + CV_PROP_RW int interpolation = INTER_LINEAR; +}; + +/** @brief Finds the geometric transform (warp) between two images in terms of the ECC criterion @cite EP08. Uses pyramids. + +@param reference Single channel reference image; CV_8U, CV_16U, CV_32F, CV_64F type. +@param sample sample image which should be warped with the final warpMatrix in +order to provide an image similar to reference, same type as reference. +@param warpMatrix floating-point \f$2\times 3\f$ or \f$3\times 3\f$ mapping matrix (warp). +@param eccParams List of the algorithm parameters. See ECCParameters for details. +@param referenceMask An optional single channel mask to indicate valid values of reference. +@param sampleMask An optional single channel mask to indicate valid values of sample. + +The function estimates the optimum transformation (warpMatrix) with respect to ECC criterion +(@cite EP08), that is + +\f[\texttt{warpMatrix} = \arg\max_{W} \texttt{ECC}(\texttt{templateImage}(x,y),\texttt{inputImage}(x',y'))\f] + +where + +\f[\begin{bmatrix} x' \\ y' \end{bmatrix} = W \cdot \begin{bmatrix} x \\ y \\ 1 \end{bmatrix}\f] + +(the equation holds with homogeneous coordinates for homography). It returns the final enhanced +correlation coefficient, that is the correlation coefficient between the template image and the +final warped input image. When a \f$3\times 3\f$ matrix is given with motionType =0, 1 or 2, the third +row is ignored. + +Unlike findHomography and estimateRigidTransform, the function findTransformECCMultiScale implements +an area-based alignment that builds on intensity similarities. In essence, the function updates the +initial transformation that roughly aligns the images. If this information is missing, the identity +warp (unity matrix) is used as an initialization. Note that if images undergo strong +displacements/rotations, an initial transformation that roughly aligns the images is necessary +(e.g., a simple euclidean/similarity transform that allows for the images showing the same image +content approximately). Use inverse warping in the second image to take an image close to the first +one, i.e. use the flag WARP_INVERSE_MAP with warpAffine or warpPerspective. See also the OpenCV +sample image_alignment.cpp that demonstrates the use of the function. Note that the function throws +an exception if algorithm does not converges. +Unlike findTransformECC, the findTransformECCMultiScale uses pyramids, making function more stable +and able to handle correctly more sophisticated cases. + +@sa +computeECC, estimateAffine2D, estimateAffinePartial2D, findHomography +*/ + +CV_EXPORTS_W double findTransformECCMultiScale(InputArray reference, + InputArray sample, + InputOutputArray warpMatrix, + const ECCParameters& eccParams = ECCParameters(), + InputArray referenceMask = noArray(), + InputArray sampleMask = noArray()); + /** @example samples/cpp/kalman.cpp An example using the standard Kalman filter */ diff --git a/modules/video/perf/perf_ecc.cpp b/modules/video/perf/perf_ecc.cpp index f81f64d8a6..87f1fbf32f 100644 --- a/modules/video/perf/perf_ecc.cpp +++ b/modules/video/perf/perf_ecc.cpp @@ -5,9 +5,14 @@ using namespace perf; CV_ENUM(MotionType, MOTION_TRANSLATION, MOTION_EUCLIDEAN, MOTION_AFFINE, MOTION_HOMOGRAPHY) CV_ENUM(ReadFlag, IMREAD_GRAYSCALE, IMREAD_COLOR) +CV_ENUM(MultiScaleFlag, false, true) typedef std::tuple TestParams; +typedef std::tuple TestParamsMS; typedef perf::TestBaseWithParam ECCPerfTest; +typedef perf::TestBaseWithParam ECCPerfTestMS; + +typedef std::tuple TestParams; PERF_TEST_P(ECCPerfTest, findTransformECC, testing::Combine(testing::Values(MOTION_TRANSLATION, MOTION_EUCLIDEAN, MOTION_AFFINE, MOTION_HOMOGRAPHY), @@ -67,4 +72,59 @@ PERF_TEST_P(ECCPerfTest, findTransformECC, } } +PERF_TEST_P(ECCPerfTestMS, findTransformECCMultiScale, + testing::Combine(testing::Values(MOTION_TRANSLATION, MOTION_EUCLIDEAN, MOTION_AFFINE, MOTION_HOMOGRAPHY), + testing::Values(false, true))) { + int transform_type = get<0>(GetParam()); + bool multiscaleFlag = get<1>(GetParam()); + + Mat img = imread(getDataPath("cv/shared/3MP.png"), IMREAD_GRAYSCALE); + Mat templateImage; + + Mat warpMat; + Mat warpGround; + + double angle; + switch (transform_type) { + case MOTION_TRANSLATION: + warpGround = (Mat_(2, 3) << 1.f, 0.f, 7.234f, 0.f, 1.f, 11.839f); + warpAffine(img, templateImage, warpGround, img.size(), INTER_LINEAR + WARP_INVERSE_MAP); + break; + case MOTION_EUCLIDEAN: + angle = CV_PI / 30; + + warpGround = (Mat_(2, 3) << (float)cos(angle), (float)-sin(angle), 12.123f, (float)sin(angle), + (float)cos(angle), 14.789f); + warpAffine(img, templateImage, warpGround, img.size(), INTER_LINEAR + WARP_INVERSE_MAP); + break; + case MOTION_AFFINE: + warpGround = (Mat_(2, 3) << 0.98f, 0.03f, 15.523f, -0.02f, 0.95f, 10.456f); + warpAffine(img, templateImage, warpGround, img.size(), INTER_LINEAR + WARP_INVERSE_MAP); + break; + case MOTION_HOMOGRAPHY: + warpGround = (Mat_(3, 3) << 0.98f, 0.03f, 15.523f, -0.02f, 0.95f, 10.456f, 0.0002f, 0.0003f, 1.f); + warpPerspective(img, templateImage, warpGround, img.size(), INTER_LINEAR + WARP_INVERSE_MAP); + break; + } + + TEST_CYCLE() { + if (transform_type < 3) + warpMat = Mat::eye(2, 3, CV_32F); + else + warpMat = Mat::eye(3, 3, CV_32F); + + if(multiscaleFlag) { + ECCParameters params; + params.criteria = cv::TermCriteria(TermCriteria::COUNT + TermCriteria::EPS, 5, -1); + params.motionType = transform_type; + params.itersPerLevel = {1, 2, 2, 2}; + findTransformECCMultiScale(templateImage, img, warpMat, params); + } + else { + findTransformECC(templateImage, img, warpMat, transform_type, + TermCriteria(TermCriteria::COUNT + TermCriteria::EPS, 5, -1)); + } + } + SANITY_CHECK_NOTHING(); +} } // namespace opencv_test diff --git a/modules/video/src/eccms.cpp b/modules/video/src/eccms.cpp new file mode 100644 index 0000000000..963763d5d3 --- /dev/null +++ b/modules/video/src/eccms.cpp @@ -0,0 +1,885 @@ +// This file is part of OpenCV project. +// It is subject to the license terms in the LICENSE file found in the top-level directory +// of this distribution and at http://opencv.org/license.html + +#include "precomp.hpp" + +/****************************************************************************************\ +* Image Alignment (ECC algorithm, pyramidal version) * +\****************************************************************************************/ + +namespace cv { +typedef std::vector MatPyramid; + +template struct MotionTraits {}; + +template<> struct MotionTraits { + enum { paramAmount = 2 }; + static inline void tailHandlerGetCoord(float& sx, + float& sy, + float& denominator, + int col, + float numeratorX0, + float numeratorY0, + float /*denominator0*/, + float /*a00*/, + float /*a10*/, + float /*a20*/) + { + denominator = 0; + sx = (numeratorX0 + col); + sy = numeratorY0; + } + template + static constexpr std::array fillJacobian(int /*col*/, int /*row*/, float/*sx*/, float/*sy*/, float fVal, + elemtype gx, elemtype gy, float /*a00*/, float /*a10*/, + float/*denominator*/) { +#define GX (fVal * gx) +#define GY (fVal * gy) + return std::array{GX, GY}; +#undef GX +#undef GY + } +}; + +template<> struct MotionTraits { + enum { paramAmount = 3 }; + static inline void tailHandlerGetCoord(float& sx, + float& sy, + float& denominator, + int col, + float numeratorX0, + float numeratorY0, + float /*denominator0*/, + float a00, + float a10, + float /*a20*/) + { + denominator = 0; + sx = (numeratorX0 + a00 * col); + sy = (numeratorY0 + a10 * col); + } + + template + static constexpr std::array fillJacobian(int col, int row, float/*sx*/, float/*sy*/, float fVal, + elemtype gx, elemtype gy, float a00, float a10, + float/*denominator*/) { +#define GX (fVal * gx) +#define GY (fVal * gy) +#define HATX (-col * a10 - row * a00) +#define HATY (col * a00 - row * a10) +#define GZ (GX * HATX + GY * HATY) + return std::array{GZ, GX, GY}; +#undef GX +#undef GY +#undef HATX +#undef HATY +#undef GZ + } +}; + +template<> struct MotionTraits { + enum { paramAmount = 6}; + static inline void tailHandlerGetCoord(float& sx, + float& sy, + float& denominator, + int col, + float numeratorX0, + float numeratorY0, + float /*denominator0*/, + float a00, + float a10, + float /*a20*/) + { + denominator = 0; + sx = (numeratorX0 + a00 * col); + sy = (numeratorY0 + a10 * col); + } + + template + static constexpr std::array fillJacobian(int col, int row, float/*sx*/, float/*sy*/, float fVal, + elemtype gx, elemtype gy, float /*a00*/, float /*a10*/, + float/*denominator*/) { +#define GX (fVal * gx) +#define GY (fVal * gy) + return std::array{GX * col, GY * col, GX * row, GY * row, GX, GY}; +#undef GX +#undef GY + } +}; + +template<> struct MotionTraits { + enum { paramAmount = 8}; + static inline void tailHandlerGetCoord(float& sx, + float& sy, + float& denominator, + int col, + float numeratorX0, + float numeratorY0, + float denominator0, + float a00, + float a10, + float a20) + { + denominator = 1.f / (col * a20 + denominator0); + sx = (numeratorX0 + a00 * col) * denominator; + sy = (numeratorY0 + a10 * col) * denominator; + } + + template + static constexpr std::array fillJacobian(int col, int row, float sx, float sy, float fVal, + elemtype gx, elemtype gy, float/*a00*/, float/*a10*/, + float denominator) { +#define GX (fVal * float(gx) * denominator) +#define GY (fVal * float(gy) * denominator) +#define GZ (-(GX * sx + GY * sy)) + return std::array{GX * col, GY * col, GZ * col, GX * row, GY * row, GZ * row, GX, GY}; +#undef GX +#undef GY +#undef GZ + } +}; + +inline void reinterpret(Mat& mat, int newdepth) { + mat.flags = (mat.flags & ~CV_MAT_DEPTH_MASK) | newdepth; +} + +template +class constexprForClass +{ +public: + static inline void execute(F&& fVal) { + constexprForClass::execute(std::forward(fVal)); + fVal(N-1); + } +}; + +template +class constexprForClass<0, F> +{ +public: + static inline void execute(F&&) {} +}; + +template +void constexprFor(F&& fVal) { + constexprForClass::execute(std::forward(fVal)); +} +template +class constexprForUpperTriangleClassOneRow +{ +public: + static inline void execute(F&& fVal) { + constexprForUpperTriangleClassOneRow::execute(std::forward(fVal)); + fVal(R, R + C - 1); + } +}; + +template +class constexprForUpperTriangleClassOneRow +{ +public: + static inline void execute(F&&) {} +}; + +template +class constexprForUpperTriangleClass +{ +public: + static inline void execute(F&& fVal) { + constexprForUpperTriangleClass::execute(std::forward(fVal)); + constexprForUpperTriangleClassOneRow::execute(std::forward(fVal)); + } +}; + +template +class constexprForUpperTriangleClass<0, D, F> +{ +public: + static inline void execute(F&&) {} +}; + +template +void constexprForUpperTriangle(F&& fVal) { + constexprForUpperTriangleClass::execute(std::forward(fVal)); +} + +template +constexpr int hessianRowStart(int row) { + return row == 0 ? 0 : (MotionTraits::paramAmount - row + 1 + hessianRowStart(row - 1)); +} + +template +static double imageHessianProjECC(const Mat& map, + const Mat& sampleWithGrad, + const Mat& ref, + double& sampSum, + double& sampSqSum, + double& refSum, + double& refSqSum, + int& nz, + Mat& hessian, + Mat& sampleProj, + Mat& refProj, + int deltaY, + int interpolation) { + static_assert(std::is_same::value, "imageHessianProjECC: f16 is not supported yet"); +#define HESSIAN_PARAMS (MotionTraits::paramAmount) + + CV_Assert(map.type() == CV_64F); + CV_Assert(interpolation == INTER_NEAREST || interpolation == INTER_LINEAR); + CV_Assert(hessian.type() == CV_64F && sampleProj.type() == CV_64F && refProj.type() == CV_64F); + if (sampleProj.size() != Size(1, HESSIAN_PARAMS) || refProj.size() != Size(1, HESSIAN_PARAMS)) { + CV_Error(Error::BadImageSize, format("imageHessianProjECC: Wrong sample projection/reference projection size. 1x%d expected", HESSIAN_PARAMS)); + } + if (hessian.size() != Size(HESSIAN_PARAMS, HESSIAN_PARAMS)) { + CV_Error(Error::BadImageSize, format("imageHessianProjECC: Wrong hessian size. %dx%d expected", HESSIAN_PARAMS, HESSIAN_PARAMS)); + } + if (!map.isContinuous()) { + CV_Error(Error::BadStep, "imageHessianProjECC: Map should be continuous"); + } + if (std::is_same::value) { + CV_Assert(sampleWithGrad.type() == CV_32FC4 && ref.type() == CV_32FC2); + } + + int hr = ref.rows; + int wr = ref.cols; + int hs = sampleWithGrad.rows; + int ws = sampleWithGrad.cols; + unsigned int ycond = hs - (INTER_LINEAR ? 1 : 0); + unsigned int xcond = ws - (INTER_LINEAR ? 1 : 0); + + hessian = Mat::zeros(hessian.size(), hessian.type()); + sampleProj = Mat::zeros(sampleProj.size(), sampleProj.type()); + refProj = Mat::zeros(refProj.size(), refProj.type()); + + const int MAX_STRIPES = 128; + int stripesAmount = std::min(MAX_STRIPES, hr / deltaY); + std::vector::paramAmount, MotionTraits::paramAmount> > hessPs(stripesAmount); + std::vector::paramAmount> > iprojs(stripesAmount); + std::vector::paramAmount> > tprojs(stripesAmount); + std::vector::paramAmount> > projSubs(stripesAmount); + std::vector correlations(stripesAmount, 0.); + + std::vector sampSums(stripesAmount, 0); + std::vector sampSqSums(stripesAmount, 0); + std::vector refSums(stripesAmount, 0); + std::vector refSqSums(stripesAmount, 0); + std::vector nzs(stripesAmount, 0); + std::vector sampMaskedSums(stripesAmount, 0); + std::vector refMaskedSums(stripesAmount, 0); + + double a00 = map.at(0, 0); + double a01 = map.at(0, 1); + double a02 = map.at(0, 2); + double a10 = map.at(1, 0); + double a11 = map.at(1, 1); + double a12 = map.at(1, 2); + double a20 = 0; + double a21 = 0; + double a22 = 0; + if (motionType == MOTION_HOMOGRAPHY) { + a20 = map.at(2, 0); + a21 = map.at(2, 1); + a22 = map.at(2, 2); + } + + const elemtype* samplePtr0 = sampleWithGrad.ptr(0); + + parallel_for_(Range(0, stripesAmount), [&](const Range& range) { + int stripeIdx = range.start; + int ystart = (hr * stripeIdx) / stripesAmount; + ystart = roundUp(ystart, deltaY); + int yend = (hr * (range.end)) / stripesAmount; + // we don't store intermediate jacobian; instead, we iteratively update Hessian, sampleProj and refProj + for (int y = ystart; y < yend; y += deltaY) { + const elemtype* refPtr = ref.ptr(y); + + std::array hessPcache{}; + std::array iprojCache{}; + std::array tprojCache{}; + std::array projSubCache{}; + + const float numeratorX0 = y * (float)a01 + (float)a02; + const float numeratorY0 = y * (float)a11 + (float)a12; + const float denominator0 = y * (float)a21 + (float)a22; + int x = 0; + for (; x < wr; x++) { //Tail handler + float sx, sy, denominator; + MotionTraits::tailHandlerGetCoord(sx, sy, denominator, x, numeratorX0, numeratorY0, + denominator0, (float)a00, (float)a10, (float)a20); + const unsigned int x0 = (interpolation == INTER_LINEAR) ? static_cast(std::floor(sx)) : saturate_cast(sx); + const unsigned int y0 = (interpolation == INTER_LINEAR) ? static_cast(std::floor(sy)) : saturate_cast(sy); + if(interpolation == INTER_LINEAR && (static_cast(x0 < xcond) & static_cast(y0 < ycond)) == 0) + continue; + if (interpolation == INTER_NEAREST && (static_cast(x0 < xcond) & static_cast(y0 < ycond)) == 0) + continue; + float sampleVal = 0; + float gx = 0; + float gy = 0; + float fVal = 0; + if(interpolation == INTER_LINEAR) { + const int x1 = x0 + 1; + const int y1 = y0 + 1; + + const float dx = sx - x0; + const float dy = sy - y0; + + const float p00_val = samplePtr0[4 * y0 * ws + 4 * x0]; + const float p01_val = samplePtr0[4 * y0 * ws + 4 * x1]; + const float p10_val = samplePtr0[4 * y1 * ws + 4 * x0]; + const float p11_val = samplePtr0[4 * y1 * ws + 4 * x1]; + const float p0_val = p00_val * (1.0f - dx) + p01_val * dx; + const float p1_val = p10_val * (1.0f - dx) + p11_val * dx; + + const float p00_gx = samplePtr0[4 * y0 * ws + 4 * x0 + 1]; + const float p01_gx = samplePtr0[4 * y0 * ws + 4 * x1 + 1]; + const float p10_gx = samplePtr0[4 * y1 * ws + 4 * x0 + 1]; + const float p11_gx = samplePtr0[4 * y1 * ws + 4 * x1 + 1]; + const float p0_gx = p00_gx * (1.0f - dx) + p01_gx * dx; + const float p1_gx = p10_gx * (1.0f - dx) + p11_gx * dx; + + const float p00_gy = samplePtr0[4 * y0 * ws + 4 * x0 + 2]; + const float p01_gy = samplePtr0[4 * y0 * ws + 4 * x1 + 2]; + const float p10_gy = samplePtr0[4 * y1 * ws + 4 * x0 + 2]; + const float p11_gy = samplePtr0[4 * y1 * ws + 4 * x1 + 2]; + const float p0_gy = p00_gy * (1.0f - dx) + p01_gy * dx; + const float p1_gy = p10_gy * (1.0f - dx) + p11_gy * dx; + + const float p00_mask = samplePtr0[4 * y0 * ws + 4 * x0 + 2] == 0.f ? 0.f : 1.f; + const float p01_mask = samplePtr0[4 * y0 * ws + 4 * x1 + 2] == 0.f ? 0.f : 1.f; + const float p10_mask = samplePtr0[4 * y1 * ws + 4 * x0 + 2] == 0.f ? 0.f : 1.f; + const float p11_mask = samplePtr0[4 * y1 * ws + 4 * x1 + 2] == 0.f ? 0.f : 1.f; + + sampleVal = p0_val * (1.0f - dy) + p1_val * dy; + gx = p0_gx * (1.0f - dy) + p1_gx * dy; + gy = p0_gy * (1.0f - dy) + p1_gy * dy; + fVal = p00_mask * p01_mask * p10_mask * p11_mask; + } + else { // if(interpolation == INTER_NEAREST) + const elemtype* samplePtr = samplePtr0 + y0 * (ws * 4) + x0 * 4; + sampleVal = samplePtr[0]; + gx = samplePtr[1]; + gy = samplePtr[2]; + fVal = float(samplePtr[3]) == 0.f ? 0.f : 1.f; + } + + float refVal = refPtr[2 * x]; + fVal *= float(refPtr[2 * x + 1]) == 0.f ? 0.f : 1.f; + sampleVal *= fVal; + refVal *= fVal; + sampSums[stripeIdx] += sampleVal; + sampSqSums[stripeIdx] += sampleVal * sampleVal; + refSums[stripeIdx] += refVal; + refSqSums[stripeIdx] += refVal * refVal; + nzs[stripeIdx] += (int)fVal; + sampMaskedSums[stripeIdx] += sampleVal; + refMaskedSums[stripeIdx] += refVal; + std::array jac = MotionTraits::fillJacobian(x, y, sx, sy, + fVal, gx, + gy, (float)a00, + (float)a10, denominator); + constexprForUpperTriangle([&](int row_i, int col_i) { + hessPcache[hessianRowStart(row_i) + (col_i - row_i)] += jac[row_i] * jac[col_i]; + }); + constexprFor([&](int elem) { + iprojCache[elem] += jac[elem] * sampleVal; + tprojCache[elem] += jac[elem] * refVal; + projSubCache[elem] += jac[elem] * fVal; + }); + correlations[stripeIdx] += sampleVal * refVal; + } + + constexprForUpperTriangle([&](int row, int col) { + hessPs[stripeIdx](row, col) += hessPcache[hessianRowStart(row) + (col - row)]; + }); + constexprFor([&](int elem) { + iprojs[stripeIdx][elem] += iprojCache[elem]; + tprojs[stripeIdx][elem] += tprojCache[elem]; + projSubs[stripeIdx][elem] += projSubCache[elem]; + }); + } + }); + double sampMaskedSum = 0; + double refMaskedSum = 0; + double correlation = 0; + sampSum = sampSqSum = refSum = refSqSum = nz = 0; + + for (int stripeIdx = 0; stripeIdx < stripesAmount; stripeIdx++) { + correlation += correlations[stripeIdx]; + sampSum += sampSums[stripeIdx]; + sampSqSum += sampSqSums[stripeIdx]; + refSum += refSums[stripeIdx]; + refSqSum += refSqSums[stripeIdx]; + sampMaskedSum += sampMaskedSums[stripeIdx]; + refMaskedSum += refMaskedSums[stripeIdx]; + nz += nzs[stripeIdx]; + } + double scale = nz == 0 ? 0. : 1. / nz; + double sampMean = sampSum * scale; + double refMean = refSum * scale; + correlation += nz * sampMean * refMean - sampMaskedSum * refMean - refMaskedSum * sampMean; + double* hessPtr = hessian.ptr(0); + double* sampleProjPtr = sampleProj.ptr(0); + double* refProjPtr = refProj.ptr(0); + for (int stripeIdx = 0; stripeIdx < stripesAmount; stripeIdx++) { + for (int hessNum = 0; hessNum < HESSIAN_PARAMS * HESSIAN_PARAMS; hessNum++) { + hessPtr[hessNum] += hessPs[stripeIdx].val[hessNum]; + } + for (int projNum = 0; projNum < HESSIAN_PARAMS; projNum++) { + sampleProjPtr[projNum] += iprojs[stripeIdx][projNum] - projSubs[stripeIdx][projNum] * sampMean; + refProjPtr[projNum] += tprojs[stripeIdx][projNum] - projSubs[stripeIdx][projNum] * refMean; + } + } + constexprForUpperTriangle([&](int row, int col) { + hessPtr[col * HESSIAN_PARAMS + row] = hessPtr[row * HESSIAN_PARAMS + col]; + }); + return correlation; +#undef HESSIAN_PARAMS +} + +static void updateWarpingMatrixECC(Mat& map_matrix, const Mat& update, const int motionType) { + CV_Assert(map_matrix.type() == CV_64FC1); + CV_Assert(update.type() == CV_64FC1); + + CV_Assert(motionType == MOTION_TRANSLATION || motionType == MOTION_EUCLIDEAN || motionType == MOTION_AFFINE || + motionType == MOTION_HOMOGRAPHY); + + if (motionType == MOTION_HOMOGRAPHY) + CV_Assert(map_matrix.rows == 3 && update.rows == 8); + else if (motionType == MOTION_AFFINE) + CV_Assert(map_matrix.rows == 2 && update.rows == 6); + else if (motionType == MOTION_EUCLIDEAN) + CV_Assert(map_matrix.rows == 2 && update.rows == 3); + else + CV_Assert(map_matrix.rows == 2 && update.rows == 2); + + CV_Assert(update.cols == 1); + + CV_Assert(map_matrix.isContinuous()); + CV_Assert(update.isContinuous()); + + double* mapPtr = map_matrix.ptr(0); + const double* updatePtr = update.ptr(0); + + if (motionType == MOTION_TRANSLATION) { + mapPtr[2] += updatePtr[0]; + mapPtr[5] += updatePtr[1]; + } + if (motionType == MOTION_AFFINE) { + mapPtr[0] += updatePtr[0]; + mapPtr[3] += updatePtr[1]; + mapPtr[1] += updatePtr[2]; + mapPtr[4] += updatePtr[3]; + mapPtr[2] += updatePtr[4]; + mapPtr[5] += updatePtr[5]; + } + if (motionType == MOTION_HOMOGRAPHY) { + mapPtr[0] += updatePtr[0]; + mapPtr[3] += updatePtr[1]; + mapPtr[6] += updatePtr[2]; + mapPtr[1] += updatePtr[3]; + mapPtr[4] += updatePtr[4]; + mapPtr[7] += updatePtr[5]; + mapPtr[2] += updatePtr[6]; + mapPtr[5] += updatePtr[7]; + } + if (motionType == MOTION_EUCLIDEAN) { + double new_theta = updatePtr[0]; + new_theta += asin(mapPtr[3]); + + mapPtr[2] += updatePtr[1]; + mapPtr[5] += updatePtr[2]; + mapPtr[0] = mapPtr[4] = cos(new_theta); + mapPtr[3] = sin(new_theta); + mapPtr[1] = -mapPtr[3]; + } +} + +static void optimizeECC(Mat& sampleWithGrad, + const Mat& reference, + Mat& map, + int motionType, + double* rho, + double* lastRho, + int deltaY, + int nparams, + int interpolation) { + CV_Assert(interpolation == INTER_NEAREST || interpolation == INTER_LINEAR); + + // warp-back portion of the inputImage and gradients to the coordinate space of the referenceFloat + double correlation = 0; + + // matrices needed for solving linear equation system for maximizing ECC + Mat hessian = Mat(nparams, nparams, CV_64F); + Mat hessianInv = Mat(nparams, nparams, CV_64F); + Mat sampleProjection = Mat(nparams, 1, CV_64F); + Mat referenceProjection = Mat(nparams, 1, CV_64F); + Mat sampleProjectionHessian = Mat(nparams, 1, CV_64F); + Mat errorProjection = Mat(nparams, 1, CV_64F); + Mat deltaP = Mat(nparams, 1, CV_64F); + + double sampSum; + double sampSqSum; + double referenceSum; + double referenceSqSum; + int nz; + + { // if(imageWithGrad.type() == CV_32FC4) + if (motionType == MOTION_TRANSLATION) { + correlation = imageHessianProjECC(map, + sampleWithGrad, + reference, + sampSum, + sampSqSum, + referenceSum, + referenceSqSum, + nz, + hessian, + sampleProjection, + referenceProjection, + deltaY, + interpolation); + } else if (motionType == MOTION_EUCLIDEAN) { + correlation = imageHessianProjECC(map, + sampleWithGrad, + reference, + sampSum, + sampSqSum, + referenceSum, + referenceSqSum, + nz, + hessian, + sampleProjection, + referenceProjection, + deltaY, + interpolation); + } else if (motionType == MOTION_AFFINE) { + correlation = imageHessianProjECC(map, + sampleWithGrad, + reference, + sampSum, + sampSqSum, + referenceSum, + referenceSqSum, + nz, + hessian, + sampleProjection, + referenceProjection, + deltaY, + interpolation); + } else { + correlation = imageHessianProjECC(map, + sampleWithGrad, + reference, + sampSum, + sampSqSum, + referenceSum, + referenceSqSum, + nz, + hessian, + sampleProjection, + referenceProjection, + deltaY, + interpolation); + } + } + double scale = nz == 0 ? 0. : 1. / nz; + double sampMean = sampSum * scale; + double refMean = referenceSum * scale; + double sampStd = std::sqrt(std::max(sampSqSum * scale - sampMean * sampMean, 0.)); + double refStd = std::sqrt(std::max(referenceSqSum * scale - refMean * refMean, 0.)); + + // inverse of Hessian + hessianInv = hessian.inv(); + // calculate enhanced correlation coefficient (ECC)->rho + *lastRho = *rho; + double refNorm = std::sqrt(nz * refStd * refStd); + double sampNorm = std::sqrt(nz * sampStd * sampStd); + + *rho = correlation / (sampNorm * refNorm); + if ((bool)cvIsNaN(*rho)) { + CV_Error(Error::StsNoConv, "NaN encountered."); + } + + // calculate the parameter lambda to account for illumination variation + sampleProjectionHessian = hessianInv * sampleProjection; + const double lambdaN = (sampNorm * sampNorm) - sampleProjection.dot(sampleProjectionHessian); + const double lambdaD = correlation - referenceProjection.dot(sampleProjectionHessian); + + if (lambdaD <= 0.0) { + CV_Error(Error::StsNoConv, "The algorithm stopped before its convergence. The correlation is going to be minimized. " + "Images may be uncorrelated or non-overlapped"); + } + const double lambda = (lambdaN / lambdaD); + + // estimate the update step delta_p + errorProjection = lambda * referenceProjection - sampleProjection; + gemm(hessianInv, errorProjection, 1., noArray(), 0., deltaP); + + // update warping matrix + updateWarpingMatrixECC(map, deltaP, motionType); +} + +static Mat prepareGradients(const Mat& sample) { + CV_Assert(sample.type() == CV_32FC2 || sample.type() == CV_16FC2); + + const int ws = sample.cols; + const int hs = sample.rows; + + Mat sampleWithGrad; + int ntasks = std::min(4, hs); + + { + sampleWithGrad = Mat(hs, ws, CV_32FC4); + float* dstPtr = sampleWithGrad.ptr(); + parallel_for_(Range(0, ntasks), [&](const Range& range) { + int rowstart = range.start * hs / ntasks; + int rowend = range.end * hs / ntasks; + for (int row = rowstart; row < rowend; row++) { + const float* sampleCurLine = sample.ptr(row); + const float* samplePrevLine = sample.ptr(std::max(row - 1, 0)); + const float* sampleNextLine = sample.ptr(std::min(row + 1, hs - 1)); + float gradDivY = (row > 0 && row + 1 < hs) ? 0.5f : 0.25f; + int col = 0; + for (; col < ws; col++) { + int prevCol = std::max(col - 1, 0); + int nextCol = std::min(col + 1, ws - 1); + float gradDivX = (col > 0 && col + 1 < ws) ? 0.5f : 0.25f; + dstPtr[row * ws * 4 + col * 4] = sampleCurLine[2 * col]; + dstPtr[row * ws * 4 + col * 4 + 1] = + gradDivX * (sampleCurLine[2 * nextCol] - sampleCurLine[2 * prevCol]); + dstPtr[row * ws * 4 + col * 4 + 2] = gradDivY * (sampleNextLine[2 * col] - samplePrevLine[2 * col]); + dstPtr[row * ws * 4 + col * 4 + 3] = sampleCurLine[2 * col + 1]; + } + } + }, ntasks); + } + return sampleWithGrad; +} + +static void buildPyramidECC(InputArray inputImage, + MatPyramid& imgPyramid, + InputArray& mask, + MatPyramid& maskPyramid, + int nlevels) { + imgPyramid.resize(nlevels); + inputImage.getMat().convertTo(imgPyramid[0], CV_8UC1); + maskPyramid.resize(nlevels); + if (!mask.empty()) { + mask.getMat().convertTo(maskPyramid[0], CV_8UC1); + } + for (int pyrLevel = 0; pyrLevel < nlevels - 1; ++pyrLevel) { + Size size = Size((imgPyramid[pyrLevel].cols + 1) / 2, (imgPyramid[pyrLevel].rows + 1) / 2); + pyrDown(imgPyramid[pyrLevel], imgPyramid[pyrLevel + 1], size); + if (!mask.empty()) { + pyrDown(maskPyramid[pyrLevel], maskPyramid[pyrLevel + 1], size); + threshold(maskPyramid[pyrLevel + 1], maskPyramid[pyrLevel + 1], 254, 0xff, THRESH_BINARY); + } + } +} + +static Mat spliceWithMask(const Mat& image, const Mat& mask) { + CV_Assert(image.type() == CV_32F && (mask.empty() || mask.type() == CV_8U)); + if (!mask.empty() && image.size() != mask.size()) { + CV_Error(Error::BadImageSize, "spliceWithMask: Mask and image have to be of same size."); + } + const int hs = image.rows; + const int ws = image.cols; + + Mat result; + int ntasks = std::min(4, hs); + { + union conv_ { + uint32_t valU; + float val; + conv_() : valU(0xffffffff) {} + } conv; + result = Mat(hs, ws, CV_32FC2); + parallel_for_(Range(0, ntasks), [&](const Range& range) { + int rowstart = range.start * hs / ntasks; + int rowend = range.end * hs / ntasks; + for (int row = rowstart; row < rowend; row++) { + float* dstPtr = result.ptr(row); + const float* srcPtr = image.ptr(row); + const uint8_t* maskPtr = !mask.empty() ? mask.ptr(row) : nullptr; + int col = 0; + for (; col < ws; col++) { + dstPtr[col * 2] = srcPtr[col]; + dstPtr[col * 2 + 1] = (!maskPtr || maskPtr[col]) ? conv.val : 0; + } + } + }, ntasks); + } + return result; +} + +static void scaleWarpMatrix(Mat& warpMatrix, float scale) { + if (warpMatrix.rows == 3) { + Mat invertScaleMat = Mat(3, 3, CV_64F, 0.f); + invertScaleMat.at(0, 0) = 1.f / scale; + invertScaleMat.at(1, 1) = 1.f / scale; + invertScaleMat.at(2, 2) = 1.f; + Mat scaleMatrix = invertScaleMat.clone(); + scaleMatrix.at(0, 0) = scale; + scaleMatrix.at(1, 1) = scale; + gemm(warpMatrix, invertScaleMat, 1., noArray(), 0., warpMatrix); + gemm(scaleMatrix, warpMatrix, 1., noArray(), 0., warpMatrix); + // Normalization, internal algorithms assumes, that a22 = 1.0f + for (int mel = 0; mel < 8; mel++) { + (reinterpret_cast(warpMatrix.data))[mel] /= + (reinterpret_cast(warpMatrix.data))[8]; + } + (reinterpret_cast(warpMatrix.data))[8] = 1.f; + } else { + warpMatrix.at(0, 2) *= scale; + warpMatrix.at(1, 2) *= scale; + } +} + +static void checkParams(const MatPyramid& referencePyramid, + const MatPyramid& samplePyramid, + Mat& map, + std::vector& itersPerLevel, + const ECCParameters& eccParams) { + if (itersPerLevel.empty()) { + itersPerLevel.resize(eccParams.nlevels, eccParams.criteria.maxCount); + } + CV_Assert(eccParams.interpolation == INTER_NEAREST || eccParams.interpolation == INTER_LINEAR); + CV_Assert(static_cast(itersPerLevel.size()) == eccParams.nlevels); + for (const auto& lvl : referencePyramid) { + CV_Assert(!lvl.empty() && lvl.type() == referencePyramid[0].type()); + } + CV_Assert(!samplePyramid.empty()); + for (const auto& lvl : samplePyramid) { + CV_Assert(!lvl.empty() && lvl.type() == samplePyramid[0].type()); + } + CV_Assert(samplePyramid.size() == referencePyramid.size() && samplePyramid.size() == itersPerLevel.size()); + CV_Assert(referencePyramid.back().rows > 1 && referencePyramid.back().cols > 1 && + samplePyramid.back().rows > 1 && samplePyramid.back().cols > 1); + // If the user passed an un-initialized warpMatrix, initialize to identity + if (referencePyramid[0].type() != CV_32FC2 && referencePyramid[0].type() != CV_16FC2) { + CV_Error(Error::StsError, "Reference pyramid have to be prepared via prepareReferencePyramid function"); + } + // accept only 1-channel images + CV_Assert(samplePyramid[0].type() == CV_32FC2 || samplePyramid[0].type() != CV_16FC2); + CV_Assert(map.type() == CV_64FC1); + if (map.cols != 3 || (map.rows != 2 && map.rows != 3)) { + CV_Error(Error::BadImageSize, "warpMatrix has incorrect size"); + } + + if (eccParams.motionType != MOTION_TRANSLATION && eccParams.motionType != MOTION_EUCLIDEAN && + eccParams.motionType != MOTION_AFFINE && eccParams.motionType != MOTION_HOMOGRAPHY) { + CV_Error(Error::StsError, "Incorrect motion type"); + } + + if (eccParams.motionType == MOTION_HOMOGRAPHY && map.rows != 3) { + CV_Error(Error::BadImageSize, "warpMatrix has incorrect size"); + } + + if (!((bool)(eccParams.criteria.type & TermCriteria::COUNT) || (bool)(eccParams.criteria.type & TermCriteria::EPS))) { + CV_Error(Error::StsError, "Incorrect stop eccParams.criteria"); + } +} + +static MatPyramid prepareECCPyramid(InputArray image, + InputArray imageMask, // Can be empty + int gaussFiltSize, + int nlevels) { + MatPyramid imagePyramid, maskPyramid; + buildPyramidECC(image, imagePyramid, imageMask, maskPyramid, nlevels); + for (int lvl = 0; lvl < nlevels; lvl++) { + Mat imgFloat; + imagePyramid[lvl].convertTo(imgFloat, CV_32F, 1. / 255.); + if (gaussFiltSize != 0) { + GaussianBlur(imgFloat, imgFloat, Size(gaussFiltSize, gaussFiltSize), 0, 0); + } + imagePyramid[lvl] = spliceWithMask( + imgFloat, + (static_cast(maskPyramid.size()) > lvl && !maskPyramid[lvl].empty()) ? maskPyramid[lvl] : Mat()); + } + return imagePyramid; +} + +double findTransformECCMultiScale(InputArray reference, + InputArray sample, + InputOutputArray warpMatrixA, + const ECCParameters& eccParams, + InputArray referenceMask, + InputArray sampleMask) { + MatPyramid referencePyramid = prepareECCPyramid(reference, referenceMask, eccParams.gaussFiltSize, eccParams.nlevels); + MatPyramid samplePyramid = prepareECCPyramid(sample, sampleMask, eccParams.gaussFiltSize, eccParams.nlevels); + Mat& warpMatrix = warpMatrixA.getMatRef(); + std::vector itersPerLevelCopy = eccParams.itersPerLevel; + // If the user passed an un-initialized warpMatrix, initialize to identity + if (warpMatrix.empty()) + { + int rowCount = eccParams.motionType == MOTION_HOMOGRAPHY ? 3 : 2; + warpMatrix = Mat::eye(rowCount, 3, CV_64FC1); + } + int warpMatrixType = warpMatrix.type(); + if (warpMatrixType != CV_64FC1) + { + warpMatrix.convertTo(warpMatrix, CV_64FC1); + } + + checkParams(referencePyramid, + samplePyramid, + warpMatrix, + itersPerLevelCopy, + eccParams); + + int nparams = 0; + switch (eccParams.motionType) { + case MOTION_TRANSLATION: + nparams = MotionTraits::paramAmount; + break; + case MOTION_EUCLIDEAN: + nparams = MotionTraits::paramAmount; + break; + case MOTION_AFFINE: + nparams = MotionTraits::paramAmount; + break; + case MOTION_HOMOGRAPHY: + nparams = MotionTraits::paramAmount; + break; + default: + CV_Error(Error::StsBadArg, "Incorrect motion type"); + } + + const std::vector numberOfIterations = ((eccParams.criteria.type & TermCriteria::COUNT) != 0) + ? itersPerLevelCopy + : std::vector(eccParams.nlevels, 200); + const double terminationEPS = (bool)(eccParams.criteria.type & TermCriteria::EPS) ? eccParams.criteria.epsilon : -1; + + // Scale warp matrix multiple times to lower pyramid level + for (int pyrLevel = 0; pyrLevel < eccParams.nlevels - 1; pyrLevel++) { + scaleWarpMatrix(warpMatrix, 0.5); + } + double rho = -1; + for (int pyrLevel = eccParams.nlevels - 1; pyrLevel >= 0; --pyrLevel) { + const int hr = referencePyramid[pyrLevel].rows; + + Mat sampleWithGrad = prepareGradients(samplePyramid[pyrLevel]); + + const int LOW_SIZE = 200; + int deltaY = hr < LOW_SIZE ? 1 : 2; + + // iteratively update mapMatrix + double lastRho = -terminationEPS; + for (int i = 1; (i <= numberOfIterations[pyrLevel]) && (fabs(rho - lastRho) >= terminationEPS); i++) { + optimizeECC(sampleWithGrad, referencePyramid[pyrLevel], warpMatrix, eccParams.motionType, &rho, &lastRho, deltaY, nparams, eccParams.interpolation); + } + if (pyrLevel > 0) { + scaleWarpMatrix(warpMatrix, 2); + } + } + if(warpMatrixType != CV_64FC1) + { + warpMatrix.convertTo(warpMatrix, warpMatrixType); + } + // return final correlation coefficient + return rho; +} +}; +/* End of file. */ \ No newline at end of file diff --git a/modules/video/test/test_ecc.cpp b/modules/video/test/test_ecc.cpp index 0bb3045dbc..ccbaa65230 100644 --- a/modules/video/test/test_ecc.cpp +++ b/modules/video/test/test_ecc.cpp @@ -45,16 +45,30 @@ namespace opencv_test { namespace { -class CV_ECC_BaseTest : public cvtest::BaseTest { +PARAM_TEST_CASE(Video_ECC, int, bool) +{ + int motionType; + bool usePyramids; + virtual void SetUp() + { + motionType = GET_PARAM(0); + usePyramids = GET_PARAM(1); + } +}; + +class CV_ECC_Test : public cvtest::BaseTest { public: - CV_ECC_BaseTest(); - virtual ~CV_ECC_BaseTest(); + CV_ECC_Test(int motionType, bool usePyramids); + virtual ~CV_ECC_Test(); protected: + int motionType; + double MAX_RMS; // upper bound for RMS error + double computeRMS(const Mat& mat1, const Mat& mat2); bool isMapCorrect(const Mat& mat); - virtual bool test(const Mat) { return true; }; // single test + virtual bool test(const Mat img); bool testAllTypes(const Mat img); // run test for all supported data types (U8, U16, F32, F64) bool testAllChNum(const Mat img); // run test for all supported channels count (gray, RGB) @@ -62,25 +76,28 @@ class CV_ECC_BaseTest : public cvtest::BaseTest { bool checkMap(const Mat& map, const Mat& ground); - double MAX_RMS_ECC; // upper bound for RMS error int ntests; // number of tests per motion type int ECC_iterations; // number of iterations for ECC double ECC_epsilon; // we choose a negative value, so that // ECC_iterations are always executed TermCriteria criteria; + bool usePyramids; // use version of findTransformECC with pyramids }; -CV_ECC_BaseTest::CV_ECC_BaseTest() { - MAX_RMS_ECC = 0.1; - ntests = 3; - ECC_iterations = 50; - ECC_epsilon = -1; //-> negative value means that ECC_Iterations will be executed - criteria = TermCriteria(TermCriteria::COUNT + TermCriteria::EPS, ECC_iterations, ECC_epsilon); -} -CV_ECC_BaseTest::~CV_ECC_BaseTest() {} +CV_ECC_Test::CV_ECC_Test(int a_motionType, bool a_usePyramids) : motionType(a_motionType) + , MAX_RMS(0.1) + , ntests(3) + , ECC_iterations(50) + , ECC_epsilon(-1) + , criteria(TermCriteria::COUNT + TermCriteria::EPS, ECC_iterations, ECC_epsilon) + , usePyramids(a_usePyramids) + {} -bool CV_ECC_BaseTest::isMapCorrect(const Mat& map) { + +CV_ECC_Test::~CV_ECC_Test() {} + +bool CV_ECC_Test::isMapCorrect(const Mat& map) { bool tr = true; float mapVal; for (int i = 0; i < map.rows; i++) @@ -92,7 +109,7 @@ bool CV_ECC_BaseTest::isMapCorrect(const Mat& map) { return tr; } -double CV_ECC_BaseTest::computeRMS(const Mat& mat1, const Mat& mat2) { +double CV_ECC_Test::computeRMS(const Mat& mat1, const Mat& mat2) { CV_Assert(mat1.rows == mat2.rows); CV_Assert(mat1.cols == mat2.cols); @@ -102,13 +119,13 @@ double CV_ECC_BaseTest::computeRMS(const Mat& mat1, const Mat& mat2) { return sqrt(errorMat.dot(errorMat) / (mat1.rows * mat1.cols * mat1.channels())); } -bool CV_ECC_BaseTest::checkMap(const Mat& map, const Mat& ground) { +bool CV_ECC_Test::checkMap(const Mat& map, const Mat& ground) { if (!isMapCorrect(map)) { ts->set_failed_test_info(cvtest::TS::FAIL_INVALID_OUTPUT); return false; } - if (computeRMS(map, ground) > MAX_RMS_ECC) { + if (computeRMS(map, ground) > MAX_RMS) { ts->set_failed_test_info(cvtest::TS::FAIL_BAD_ACCURACY); ts->printf(ts->LOG, "RMS = %f", computeRMS(map, ground)); return false; @@ -116,7 +133,77 @@ bool CV_ECC_BaseTest::checkMap(const Mat& map, const Mat& ground) { return true; } -bool CV_ECC_BaseTest::testAllTypes(const Mat img) { +bool CV_ECC_Test::test(const Mat img) +{ + cv::RNG rng = ts->get_rng(); + + int progress = 0; + + for (int k = 0; k < ntests; k++) { + ts->update_context(this, k, true); + progress = update_progress(progress, k, ntests, 0); + + Mat groundMap; + switch(motionType) + { + case MOTION_TRANSLATION: + groundMap = (Mat_(2, 3) << 1, 0, (rng.uniform(10.f, 20.f)), 0, 1, (rng.uniform(10.f, 20.f))); + break; + case MOTION_EUCLIDEAN: + { + double angle = CV_PI / 30 + CV_PI * rng.uniform((double)-2.f, (double)2.f) / 180; + groundMap = (Mat_(2, 3) << cos(angle), -sin(angle), (rng.uniform(10.f, 20.f)), sin(angle), + cos(angle), (rng.uniform(10.f, 20.f))); + break; + } + case MOTION_AFFINE: + groundMap = (Mat_(2, 3) << (1 - rng.uniform(-0.05f, 0.05f)), (rng.uniform(-0.03f, 0.03f)), + (rng.uniform(10.f, 20.f)), (rng.uniform(-0.03f, 0.03f)), (1 - rng.uniform(-0.05f, 0.05f)), + (rng.uniform(10.f, 20.f))); + break; + case MOTION_HOMOGRAPHY: + groundMap = + (Mat_(3, 3) << (1 - rng.uniform(-0.05f, 0.05f)), (rng.uniform(-0.03f, 0.03f)), + (rng.uniform(10.f, 20.f)), (rng.uniform(-0.03f, 0.03f)), (1 - rng.uniform(-0.05f, 0.05f)), + (rng.uniform(10.f, 20.f)), (rng.uniform(0.0001f, 0.0003f)), (rng.uniform(0.0001f, 0.0003f)), 1.f); + break; + default: + CV_Error(Error::StsBadArg, "Incorrect motion type"); + break; + } + + Mat warpedImage; + + Mat foundMap; + if(motionType == MOTION_HOMOGRAPHY) + { + warpPerspective(img, warpedImage, groundMap, Size(200, 200), INTER_LINEAR + WARP_INVERSE_MAP); + foundMap = Mat::eye(3, 3, CV_32F); + } + else + { + warpAffine(img, warpedImage, groundMap, Size(200, 200), INTER_LINEAR + WARP_INVERSE_MAP); + foundMap = Mat((Mat_(2, 3) << 1, 0, 0, 0, 1, 0)); + } + + + if(usePyramids) + { + ECCParameters params; + params.criteria = criteria; + params.motionType = motionType; + findTransformECCMultiScale(warpedImage, img, foundMap, params); + } + else + findTransformECC(warpedImage, img, foundMap, motionType, criteria); + + if (!checkMap(foundMap, groundMap)) + return false; + } + return true; +} + +bool CV_ECC_Test::testAllTypes(const Mat img) { auto types = {CV_8U, CV_16U, CV_32F, CV_64F}; for (auto type : types) { Mat timg; @@ -127,9 +214,10 @@ bool CV_ECC_BaseTest::testAllTypes(const Mat img) { return true; } -bool CV_ECC_BaseTest::testAllChNum(const Mat img) { - if (!testAllTypes(img)) - return false; +bool CV_ECC_Test::testAllChNum(const Mat img) { + if(!usePyramids) + if (!testAllTypes(img)) + return false; Mat gray; cvtColor(img, gray, COLOR_RGB2GRAY); @@ -139,7 +227,7 @@ bool CV_ECC_BaseTest::testAllChNum(const Mat img) { return true; } -void CV_ECC_BaseTest::run(int) { +void CV_ECC_Test::run(int) { Mat img = imread(string(ts->get_data_path()) + "shared/fruits.png"); if (img.empty()) { ts->printf(ts->LOG, "test image can not be read"); @@ -155,153 +243,22 @@ void CV_ECC_BaseTest::run(int) { ts->set_failed_test_info(cvtest::TS::OK); } -class CV_ECC_Test_Translation : public CV_ECC_BaseTest { - public: - CV_ECC_Test_Translation(); - - protected: - bool test(const Mat); -}; - -CV_ECC_Test_Translation::CV_ECC_Test_Translation() {} - -bool CV_ECC_Test_Translation::test(const Mat testImg) { - cv::RNG rng = ts->get_rng(); - - int progress = 0; - - for (int k = 0; k < ntests; k++) { - ts->update_context(this, k, true); - progress = update_progress(progress, k, ntests, 0); - - Mat translationGround = (Mat_(2, 3) << 1, 0, (rng.uniform(10.f, 20.f)), 0, 1, (rng.uniform(10.f, 20.f))); - - Mat warpedImage; - - warpAffine(testImg, warpedImage, translationGround, Size(200, 200), INTER_LINEAR + WARP_INVERSE_MAP); - - Mat mapTranslation = (Mat_(2, 3) << 1, 0, 0, 0, 1, 0); - - findTransformECC(warpedImage, testImg, mapTranslation, 0, criteria); - - if (!checkMap(mapTranslation, translationGround)) - return false; - } - return true; +TEST_P(Video_ECC, accuracy) { + CV_ECC_Test test(motionType, usePyramids); + test.safe_run(); } -class CV_ECC_Test_Euclidean : public CV_ECC_BaseTest { - public: - CV_ECC_Test_Euclidean(); +INSTANTIATE_TEST_CASE_P(ECCfixtures, Video_ECC, + testing::Values(testing::make_tuple(MOTION_TRANSLATION, false), + testing::make_tuple(MOTION_TRANSLATION, true), + testing::make_tuple(MOTION_EUCLIDEAN, false), + testing::make_tuple(MOTION_EUCLIDEAN, true), + testing::make_tuple(MOTION_AFFINE, false), + testing::make_tuple(MOTION_AFFINE, true), + testing::make_tuple(MOTION_HOMOGRAPHY, false), + testing::make_tuple(MOTION_HOMOGRAPHY, true))); - protected: - bool test(const Mat); -}; - -CV_ECC_Test_Euclidean::CV_ECC_Test_Euclidean() {} - -bool CV_ECC_Test_Euclidean::test(const Mat testImg) { - cv::RNG rng = ts->get_rng(); - - int progress = 0; - for (int k = 0; k < ntests; k++) { - ts->update_context(this, k, true); - progress = update_progress(progress, k, ntests, 0); - - double angle = CV_PI / 30 + CV_PI * rng.uniform((double)-2.f, (double)2.f) / 180; - - Mat euclideanGround = (Mat_(2, 3) << cos(angle), -sin(angle), (rng.uniform(10.f, 20.f)), sin(angle), - cos(angle), (rng.uniform(10.f, 20.f))); - - Mat warpedImage; - - warpAffine(testImg, warpedImage, euclideanGround, Size(200, 200), INTER_LINEAR + WARP_INVERSE_MAP); - - Mat mapEuclidean = (Mat_(2, 3) << 1, 0, 0, 0, 1, 0); - - findTransformECC(warpedImage, testImg, mapEuclidean, 1, criteria); - - if (!checkMap(mapEuclidean, euclideanGround)) - return false; - } - return true; -} - -class CV_ECC_Test_Affine : public CV_ECC_BaseTest { - public: - CV_ECC_Test_Affine(); - - protected: - bool test(const Mat img); -}; - -CV_ECC_Test_Affine::CV_ECC_Test_Affine() {} - -bool CV_ECC_Test_Affine::test(const Mat testImg) { - cv::RNG rng = ts->get_rng(); - - int progress = 0; - for (int k = 0; k < ntests; k++) { - ts->update_context(this, k, true); - progress = update_progress(progress, k, ntests, 0); - - Mat affineGround = (Mat_(2, 3) << (1 - rng.uniform(-0.05f, 0.05f)), (rng.uniform(-0.03f, 0.03f)), - (rng.uniform(10.f, 20.f)), (rng.uniform(-0.03f, 0.03f)), (1 - rng.uniform(-0.05f, 0.05f)), - (rng.uniform(10.f, 20.f))); - - Mat warpedImage; - - warpAffine(testImg, warpedImage, affineGround, Size(200, 200), INTER_LINEAR + WARP_INVERSE_MAP); - - Mat mapAffine = (Mat_(2, 3) << 1, 0, 0, 0, 1, 0); - - findTransformECC(warpedImage, testImg, mapAffine, 2, criteria); - - if (!checkMap(mapAffine, affineGround)) - return false; - } - - return true; -} - -class CV_ECC_Test_Homography : public CV_ECC_BaseTest { - public: - CV_ECC_Test_Homography(); - - protected: - bool test(const Mat testImg); -}; - -CV_ECC_Test_Homography::CV_ECC_Test_Homography() {} - -bool CV_ECC_Test_Homography::test(const Mat testImg) { - cv::RNG rng = ts->get_rng(); - - int progress = 0; - for (int k = 0; k < ntests; k++) { - ts->update_context(this, k, true); - progress = update_progress(progress, k, ntests, 0); - - Mat homoGround = - (Mat_(3, 3) << (1 - rng.uniform(-0.05f, 0.05f)), (rng.uniform(-0.03f, 0.03f)), - (rng.uniform(10.f, 20.f)), (rng.uniform(-0.03f, 0.03f)), (1 - rng.uniform(-0.05f, 0.05f)), - (rng.uniform(10.f, 20.f)), (rng.uniform(0.0001f, 0.0003f)), (rng.uniform(0.0001f, 0.0003f)), 1.f); - - Mat warpedImage; - - warpPerspective(testImg, warpedImage, homoGround, Size(200, 200), INTER_LINEAR + WARP_INVERSE_MAP); - - Mat mapHomography = Mat::eye(3, 3, CV_32F); - - findTransformECC(warpedImage, testImg, mapHomography, 3, criteria); - - if (!checkMap(mapHomography, homoGround)) - return false; - } - return true; -} - -class CV_ECC_Test_Mask : public CV_ECC_BaseTest { +class CV_ECC_Test_Mask : public CV_ECC_Test { public: CV_ECC_Test_Mask(); @@ -309,7 +266,7 @@ class CV_ECC_Test_Mask : public CV_ECC_BaseTest { bool test(const Mat); }; -CV_ECC_Test_Mask::CV_ECC_Test_Mask() {} +CV_ECC_Test_Mask::CV_ECC_Test_Mask():CV_ECC_Test(MOTION_TRANSLATION, false) {} bool CV_ECC_Test_Mask::test(const Mat testImg) { cv::RNG rng = ts->get_rng(); @@ -368,6 +325,58 @@ bool CV_ECC_Test_Mask::test(const Mat testImg) { return true; } +class CV_ECC_BigPictureTest : public CV_ECC_Test { + public: + CV_ECC_BigPictureTest(bool a_maskedVersion) : CV_ECC_Test(MOTION_HOMOGRAPHY, true), maskedVersion(a_maskedVersion) {} + virtual ~CV_ECC_BigPictureTest() {} + protected: + void run(int); + bool maskedVersion; +}; + +void CV_ECC_BigPictureTest::run(int) +{ + Mat largeGray0 = imread(string(ts->get_data_path()) + "shared/halmosh0.jpg", IMREAD_GRAYSCALE); + Mat largeGray1; + Mat roiMask0; + Mat roiMask1; + Mat expectedRes; + bool readError = false; + if(maskedVersion) + { + largeGray1 = imread(string(ts->get_data_path()) + "shared/halmosh2.jpg", IMREAD_GRAYSCALE); + roiMask0 = imread(string(ts->get_data_path()) + "shared/halmosh0mask.png", IMREAD_GRAYSCALE); + roiMask1 = imread(string(ts->get_data_path()) + "shared/halmosh2mask.png", IMREAD_GRAYSCALE); + readError = largeGray0.empty() || largeGray1.empty() || roiMask0.empty() || roiMask1.empty(); + expectedRes = (Mat_(3, 3) << 1.0225, 0.0606, -28.6452, -0.0475, 1.0314, 11.819, 8.21e-06, -3.65e-07, 1); + } + else + { + largeGray1 = imread(string(ts->get_data_path()) + "shared/halmosh1.jpg", IMREAD_GRAYSCALE); + readError = largeGray0.empty() || largeGray1.empty(); + expectedRes = (Mat_(3, 3) << 0.9756, -0.0319, 24.685, 0.013, 0.9808, 7.7453, -2.35e-05, -9.12e-06, 1); + } + + if(readError) + { + ts->printf(ts->LOG, "test image can not be read"); + ts->set_failed_test_info(cvtest::TS::FAIL_INVALID_TEST_DATA); + return; + } + + cv::Mat found = cv::Mat::eye(3, 3, CV_32F); + constexpr int N_ITERS = 20; + constexpr double TERMINATION_EPS = 1e-6; + ECCParameters params; + params.criteria = cv::TermCriteria(cv::TermCriteria::COUNT + cv::TermCriteria::EPS, N_ITERS, TERMINATION_EPS); + params.motionType = MOTION_HOMOGRAPHY; + params.nlevels = 5; + params.itersPerLevel = {5, 10, 300, 300, 1000}; + findTransformECCMultiScale(largeGray0, largeGray1, found, params, roiMask0, roiMask1); + ASSERT_EQ(checkMap(found, expectedRes), true); + ts->set_failed_test_info(cvtest::TS::OK); +} + void testECCProperties(Mat x, float eps) { // The channels are independent Mat y = x.t(); @@ -450,26 +459,17 @@ TEST(Video_ECC_Test_Compute, bug_14657) { EXPECT_NEAR(computeECC(img, img), 1.0f, 1e-5f); } -TEST(Video_ECC_Translation, accuracy) { - CV_ECC_Test_Translation test; - test.safe_run(); -} -TEST(Video_ECC_Euclidean, accuracy) { - CV_ECC_Test_Euclidean test; - test.safe_run(); -} -TEST(Video_ECC_Affine, accuracy) { - CV_ECC_Test_Affine test; - test.safe_run(); -} -TEST(Video_ECC_Homography, accuracy) { - CV_ECC_Test_Homography test; - test.safe_run(); -} TEST(Video_ECC_Mask, accuracy) { CV_ECC_Test_Mask test; test.safe_run(); } - +TEST(Video_ECC_BigMS, accuracy) { + CV_ECC_BigPictureTest test(false); + test.safe_run(); +} +TEST(Video_ECC_BigMS_Mask, accuracy) { + CV_ECC_BigPictureTest test(true); + test.safe_run(); +} } // namespace } // namespace opencv_test