Merge pull request #29165 from vpisarev:better_bicubic

Revised bicubic interpolation kernels and updated the corresponding functions #29165

They should work more accurately and faster than the existing implementation. The PR follows the previous work (#26271 etc.) that brought faster bilinear kernels.

### Pull Request Readiness Checklist

See details at https://github.com/opencv/opencv/wiki/How_to_contribute#making-a-good-pull-request

- [x] I agree to contribute to the project under Apache 2 License.
- [x] To the best of my knowledge, the proposed patch is not based on a code under GPL or another license that is incompatible with OpenCV
- [x] The PR is proposed to the proper branch
- [x] There is a reference to the original bug report and related work
- [x] There is accuracy test, performance test and test data in opencv_extra repository, if applicable
      Patch to opencv_extra has the same branch name.
- [ ] The feature is well documented and sample code can be built with the project CMake

modules/imgproc/perf/opencl/perf_imgwarp.cpp

@@ -73,10 +73,6 @@ OCL_PERF_TEST_P(WarpAffineFixture, WarpAffine,
     const Size srcSize = get<0>(params);
     const int type = get<1>(params), interpolation = get<2>(params);
 
-    // BUG: OpenCL and CPU version diverges a bit
-    // Ticket: https://github.com/opencv/opencv/issues/26235
-    const double eps = CV_MAT_DEPTH(type) <= CV_32S ? 2 : interpolation == INTER_CUBIC ? 2e-3 : 3e-2;
-
     checkDeviceMaxMemoryAllocSize(srcSize, type);
 
     UMat src(srcSize, type), dst(srcSize, type);
@@ -84,7 +80,16 @@ OCL_PERF_TEST_P(WarpAffineFixture, WarpAffine,
 
     OCL_TEST_CYCLE() cv::warpAffine(src, dst, M, srcSize, interpolation);
 
-    SANITY_CHECK(dst, eps);
+    if (interpolation == INTER_CUBIC)
+    {
+        SANITY_CHECK_NOTHING();
+    }
+    else
+    {
+        // Ticket: https://github.com/opencv/opencv/issues/26235
+        const double eps = CV_MAT_DEPTH(type) <= CV_32S ? 2 : interpolation == INTER_CUBIC ? 2e-3 : 3e-2;
+        SANITY_CHECK(dst, eps);
+    }
 }
 
 ///////////// WarpPerspective ////////////////////////
@@ -107,7 +112,6 @@ OCL_PERF_TEST_P(WarpPerspectiveFixture, WarpPerspective,
     const WarpPerspectiveParams params = GetParam();
     const Size srcSize = get<0>(params);
     const int type = get<1>(params), interpolation = get<2>(params);
-    const double eps = CV_MAT_DEPTH(type) <= CV_32S ? 1 : 1e-4;
 
     checkDeviceMaxMemoryAllocSize(srcSize, type);
 
@@ -116,7 +120,15 @@ OCL_PERF_TEST_P(WarpPerspectiveFixture, WarpPerspective,
 
     OCL_TEST_CYCLE() cv::warpPerspective(src, dst, M, srcSize, interpolation);
 
-    SANITY_CHECK(dst, eps);
+    if (interpolation == INTER_CUBIC)
+    {
+        SANITY_CHECK_NOTHING();
+    }
+    else
+    {
+        const double eps = CV_MAT_DEPTH(type) <= CV_32S ? 1 : 1e-4;
+        SANITY_CHECK(dst, eps);
+    }
 }
 
 ///////////// Resize ////////////////////////

modules/imgproc/perf/perf_remap.cpp

@@ -64,7 +64,14 @@ PERF_TEST_P( TestRemap, Remap,
     int runs = (sz.width <= 640) ? 3 : 1;
     TEST_CYCLE_MULTIRUN(runs) remap(src, dst, map1, map2, inter_type);
 
-    SANITY_CHECK(dst);
+    if (inter_type == INTER_CUBIC)
+    {
+        SANITY_CHECK_NOTHING();
+    }
+    else
+    {
+        SANITY_CHECK(dst);
+    }
 }
 
 } // namespace

modules/imgproc/src/imgwarp.cpp

@@ -55,14 +55,162 @@
 #include "opencv2/core/softfloat.hpp"
 #include "imgwarp.hpp"
 
+#include "warp_common.hpp"
 #include "warp_kernels.simd.hpp"
 #include "warp_kernels.simd_declarations.hpp"
 
-using namespace cv;
-
 namespace cv
 {
 
+///////////////////////////// new-style kernel-base image warping without tables ////////////////////////
+
+ImgWarpFunc getImgWarpFunc(int type, int interpolation)
+{
+    if (interpolation == INTER_CUBIC) {
+        CV_CPU_DISPATCH(getBicubicWarpFunc_, (type), CV_CPU_DISPATCH_MODES_ALL);
+    }
+    return (ImgWarpFunc)nullptr;
+}
+
+static bool genericWarp(const Mat& src, const Mat& M_, const Mat& mapx, const Mat& mapy, Mat& dst,
+                        int interpolation, int borderType_, const Scalar& borderValue, bool relativeMap_)
+{
+    constexpr int MAX_CHANNELS = 4;
+    int srctype = src.type();
+
+    ImgWarpFunc func = getImgWarpFunc(srctype, interpolation);
+    if (!func) {
+        return false; // not implemented;
+    }
+
+    CV_Assert(src.type() == dst.type());
+    CV_Assert(src.channels() <= MAX_CHANNELS);
+    Matx33f Mdata(1.f, 0.f, 0.f, 0.f, 1.f, 0.f, 0.f, 0.f, 1.f);
+
+    if (!M_.empty()) {
+        if (!mapx.empty() || !mapy.empty()) {
+            CV_Error(Error::StsBadArg, "Either M or both mapx and mapy must be empty");
+        }
+        CV_Assert(M_.size() == Size(3, 2) || M_.size() == Size(3, 3));
+        CV_Assert(M_.type() == CV_32F || M_.type() == CV_64F);
+        Mat Mcopy(M_.size(), CV_32F, Mdata.val);
+        M_.convertTo(Mcopy, CV_32F);
+    } else {
+        if (mapx.empty()) {
+            CV_Error(Error::StsBadArg, "When M is empty, mapx and/or mapy must be non-empty");
+        }
+        CV_Assert(mapx.size() == dst.size());
+        if (mapy.empty()) {
+            CV_Assert(mapx.type() == CV_32FC2);
+        } else {
+            if (mapx.type() == mapy.type()) {
+                CV_Assert(mapx.type() == CV_32FC1);
+            } else {
+                CV_Assert(mapx.type() == CV_16SC2 && (mapy.type() == CV_16UC1 || mapy.type() == CV_16SC1));
+            }
+            CV_Assert(mapx.size() == mapy.size());
+        }
+    }
+
+    int64_t borderValBuf_[MAX_CHANNELS];
+    scalarToRawData(borderValue, borderValBuf_, src.type());
+
+    parallel_for_(Range(0, dst.rows), [&](const Range& range) {
+        constexpr float FIXPT_SCALE = 1.f/32;
+        constexpr int BLOCK_SIZE = 128;
+        const uint8_t* srcdata = src.data;
+        const float* fparams = nullptr;
+        size_t srcstep = src.step;
+        Size srcsize = src.size();
+        int dstcols = dst.cols;
+        int borderType = borderType_;
+        size_t bpp = src.elemSize();
+
+        bool warping = !M_.empty();
+        bool warpAffine = warping && M_.rows == 2;
+        bool warpPerspective = warping && M_.rows == 3;
+
+        bool relativeMap = relativeMap_ && !warping;
+        float relscale = relativeMap ? 1.f : 0.f;
+        bool interleavedmap = !warping && (mapx.type() == CV_32FC2);
+        bool planarmap = !warping && (mapx.type() == mapy.type());
+        bool fixedmap = !warping && (mapx.type() == CV_16SC2);
+
+        CV_Assert(warping || interleavedmap || planarmap || fixedmap);
+
+        Matx33f M = Mdata;
+        const uint8_t* borderValBuf = (const uint8_t*)borderValBuf_;
+        float M_xx = M(0, 0), M_yx = M(1, 0), M_zx = M(2, 0);
+        float xbuf[BLOCK_SIZE], ybuf[BLOCK_SIZE];
+        const float* xbufptr = xbuf;
+        const float* ybufptr = ybuf;
+
+        for (int y = range.start; y < range.end; y++) {
+            uint8_t* dstptr = dst.ptr(y);
+            float M_x = float(y)*M(0, 1) + M(0, 2);
+            float M_y = float(y)*M(1, 1) + M(1, 2);
+            float M_z = float(y)*M(2, 1) + M(2, 2);
+            const Vec2s* xysptr = fixedmap ? mapx.ptr<Vec2s>(y) : nullptr;
+            const ushort* idxptr = fixedmap ? mapy.ptr<ushort>(y) : nullptr;
+            const Vec2f* xyfptr = interleavedmap ? mapx.ptr<Vec2f>(y) : nullptr;
+            const float* xfptr = planarmap ? mapx.ptr<float>(y) : nullptr;
+            const float* yfptr = planarmap ? mapy.ptr<float>(y) : nullptr;
+            float y0 = relativeMap ? float(y) : 0.f;
+
+            for (int x = 0; x < dstcols; x += BLOCK_SIZE) {
+                int blocksize = std::min(BLOCK_SIZE, dstcols - x);
+                if (warpAffine) {
+                    for (int dx = 0; dx < blocksize; dx++) {
+                        float xf = float(x + dx);
+                        xbuf[dx] = M_x + M_xx*xf;
+                        ybuf[dx] = M_y + M_yx*xf;
+                    }
+                } else if (warpPerspective) {
+                    for (int dx = 0; dx < blocksize; dx++) {
+                        double xf = double(x + dx);
+                        double invz = 1./(M_z + M_zx*xf);
+                        xbuf[dx] = float((M_x + M_xx*xf)*invz);
+                        ybuf[dx] = float((M_y + M_yx*xf)*invz);
+                    }
+                } else if (fixedmap) {
+                    for (int dx = 0; dx < blocksize; dx++) {
+                        Vec2s xy = xysptr[x + dx];
+                        ushort idx = idxptr[x + dx];
+                        float xf = float(xy[0]) + (idx & 31)*FIXPT_SCALE + (x + dx)*relscale;
+                        float yf = float(xy[1]) + ((idx >> 5) & 31)*FIXPT_SCALE + y0;
+                        xbuf[dx] = xf;
+                        ybuf[dx] = yf;
+                    }
+                } else if (interleavedmap) {
+                    for (int dx = 0; dx < blocksize; dx++) {
+                        Vec2f xy = xyfptr[x + dx];
+                        float xf = xy[0] + (x + dx)*relscale;
+                        float yf = xy[1] + y0;
+                        xbuf[dx] = xf;
+                        ybuf[dx] = yf;
+                    }
+                } else if (relativeMap) {
+                    for (int dx = 0; dx < blocksize; dx++) {
+                        float xf = xfptr[x + dx] + (x + dx)*relscale;
+                        float yf = yfptr[x + dx] + y0;
+                        xbuf[dx] = xf;
+                        ybuf[dx] = yf;
+                    }
+                } else {
+                    // planar absolute map — just use it as-is without copying
+                    xbufptr = xfptr + x;
+                    ybufptr = yfptr + x;
+                }
+
+                func(xbufptr, ybufptr, blocksize, srcdata, srcstep, srcsize,
+                     dstptr + x*bpp, fparams, borderType, borderValBuf);
+            }
+        }
+    });
+
+    return true;
+}
+
 /************** interpolation formulas and tables ***************/
 
 const int INTER_REMAP_COEF_BITS=15;
@@ -78,9 +226,6 @@ static short BilinearTab_iC4_buf[INTER_TAB_SIZE2+2][2][8];
 static short (*BilinearTab_iC4)[2][8] = (short (*)[2][8])alignPtr(BilinearTab_iC4_buf, 16);
 #endif
 
-static float BicubicTab_f[INTER_TAB_SIZE2][4][4];
-static short BicubicTab_i[INTER_TAB_SIZE2][4][4];
-
 static float Lanczos4Tab_f[INTER_TAB_SIZE2][8][8];
 static short Lanczos4Tab_i[INTER_TAB_SIZE2][8][8];
 
@@ -90,16 +235,6 @@ static inline void interpolateLinear( float x, float* coeffs )
     coeffs[1] = x;
 }
 
-static inline void interpolateCubic( float x, float* coeffs )
-{
-    const float A = -0.75f;
-
-    coeffs[0] = ((A*(x + 1) - 5*A)*(x + 1) + 8*A)*(x + 1) - 4*A;
-    coeffs[1] = ((A + 2)*x - (A + 3))*x*x + 1;
-    coeffs[2] = ((A + 2)*(1 - x) - (A + 3))*(1 - x)*(1 - x) + 1;
-    coeffs[3] = 1.f - coeffs[0] - coeffs[1] - coeffs[2];
-}
-
 static inline void interpolateLanczos4( float x, float* coeffs )
 {
     static const double s45 = 0.70710678118654752440084436210485;
@@ -138,8 +273,7 @@ static void initInterTab1D(int method, float* tab, int tabsz)
     }
     else if( method == INTER_CUBIC )
     {
-        for( int i = 0; i < tabsz; i++, tab += 4 )
-            interpolateCubic( i*scale, tab );
+        ; // pass
     }
     else if( method == INTER_LANCZOS4 )
     {
@@ -154,13 +288,15 @@ static void initInterTab1D(int method, float* tab, int tabsz)
 static const void* initInterTab2D( int method, bool fixpt )
 {
     static bool inittab[INTER_MAX+1] = {false};
+
+    if( method == INTER_CUBIC )
+        return nullptr;
+
     float* tab = 0;
     short* itab = 0;
     int ksize = 0;
     if( method == INTER_LINEAR )
         tab = BilinearTab_f[0][0], itab = BilinearTab_i[0][0], ksize=2;
-    else if( method == INTER_CUBIC )
-        tab = BicubicTab_f[0][0], itab = BicubicTab_i[0][0], ksize=4;
     else if( method == INTER_LANCZOS4 )
         tab = Lanczos4Tab_f[0][0], itab = Lanczos4Tab_i[0][0], ksize=8;
     else
@@ -232,8 +368,6 @@ static bool initAllInterTab2D()
 {
     return  initInterTab2D( INTER_LINEAR, false ) &&
             initInterTab2D( INTER_LINEAR, true ) &&
-            initInterTab2D( INTER_CUBIC, false ) &&
-            initInterTab2D( INTER_CUBIC, true ) &&
             initInterTab2D( INTER_LANCZOS4, false ) &&
             initInterTab2D( INTER_LANCZOS4, true );
 }
@@ -845,111 +979,6 @@ static void remapBilinear( const Mat& _src, Mat& _dst, const Mat& _xy,
 }
 
 
-template<class CastOp, typename AT, int ONE, bool isRelative>
-static void remapBicubic( const Mat& _src, Mat& _dst, const Mat& _xy,
-                          const Mat& _fxy, const void* _wtab,
-                          int borderType, const Scalar& _borderValue, const Point& _offset )
-{
-    typedef typename CastOp::rtype T;
-    typedef typename CastOp::type1 WT;
-    Size ssize = _src.size(), dsize = _dst.size();
-    const int cn = _src.channels();
-    const AT* wtab = (const AT*)_wtab;
-    const T* S0 = _src.ptr<T>();
-    size_t sstep = _src.step/sizeof(S0[0]);
-    T cval[CV_CN_MAX];
-    CastOp castOp;
-
-    for(int k = 0; k < cn; k++ )
-        cval[k] = saturate_cast<T>(_borderValue[k & 3]);
-
-    int borderType1 = borderType != BORDER_TRANSPARENT ? borderType : BORDER_REFLECT_101;
-
-    unsigned width1 = std::max(ssize.width-3, 0), height1 = std::max(ssize.height-3, 0);
-
-    if( _dst.isContinuous() && _xy.isContinuous() && _fxy.isContinuous() && !isRelative )
-    {
-        dsize.width *= dsize.height;
-        dsize.height = 1;
-    }
-
-    for(int dy = 0; dy < dsize.height; dy++ )
-    {
-        T* D = _dst.ptr<T>(dy);
-        const short* XY = _xy.ptr<short>(dy);
-        const ushort* FXY = _fxy.ptr<ushort>(dy);
-        const int off_y = isRelative ? (_offset.y+dy) : 0;
-        for(int dx = 0; dx < dsize.width; dx++, D += cn )
-        {
-            const int off_x = isRelative ? (_offset.x+dx) : 0;
-            int sx = XY[dx*2]-1+off_x, sy = XY[dx*2+1]-1+off_y;
-            const AT* w = wtab + FXY[dx]*16;
-            if( (unsigned)sx < width1 && (unsigned)sy < height1 )
-            {
-                const T* S = S0 + sy*sstep + sx*cn;
-                for(int k = 0; k < cn; k++ )
-                {
-                    WT sum = S[0]*w[0] + S[cn]*w[1] + S[cn*2]*w[2] + S[cn*3]*w[3];
-                    S += sstep;
-                    sum += S[0]*w[4] + S[cn]*w[5] + S[cn*2]*w[6] + S[cn*3]*w[7];
-                    S += sstep;
-                    sum += S[0]*w[8] + S[cn]*w[9] + S[cn*2]*w[10] + S[cn*3]*w[11];
-                    S += sstep;
-                    sum += S[0]*w[12] + S[cn]*w[13] + S[cn*2]*w[14] + S[cn*3]*w[15];
-                    S -= sstep * 3 - 1;
-                    D[k] = castOp(sum);
-                }
-            }
-            else
-            {
-                int x[4], y[4];
-                if( borderType == BORDER_TRANSPARENT &&
-                    ((unsigned)(sx+1) >= (unsigned)ssize.width ||
-                    (unsigned)(sy+1) >= (unsigned)ssize.height) )
-                    continue;
-
-                if( borderType1 == BORDER_CONSTANT &&
-                    (sx >= ssize.width || sx+4 <= 0 ||
-                    sy >= ssize.height || sy+4 <= 0))
-                {
-                    for(int k = 0; k < cn; k++ )
-                        D[k] = cval[k];
-                    continue;
-                }
-
-                for(int i = 0; i < 4; i++ )
-                {
-                    x[i] = borderInterpolate(sx + i, ssize.width, borderType1)*cn;
-                    y[i] = borderInterpolate(sy + i, ssize.height, borderType1);
-                }
-
-                for(int k = 0; k < cn; k++, S0++, w -= 16 )
-                {
-                    WT cv = cval[k], sum = cv*ONE;
-                    for(int i = 0; i < 4; i++, w += 4 )
-                    {
-                        int yi = y[i];
-                        if( yi < 0 )
-                            continue;
-                        const T* S = S0 + yi*sstep;
-                        if( x[0] >= 0 )
-                            sum += (S[x[0]] - cv)*w[0];
-                        if( x[1] >= 0 )
-                            sum += (S[x[1]] - cv)*w[1];
-                        if( x[2] >= 0 )
-                            sum += (S[x[2]] - cv)*w[2];
-                        if( x[3] >= 0 )
-                            sum += (S[x[3]] - cv)*w[3];
-                    }
-                    D[k] = castOp(sum);
-                }
-                S0 -= cn;
-            }
-        }
-    }
-}
-
-
 template<class CastOp, typename AT, int ONE, bool isRelative>
 static void remapLanczos4( const Mat& _src, Mat& _dst, const Mat& _xy,
                            const Mat& _fxy, const void* _wtab,
@@ -1395,7 +1424,6 @@ void cv::remap( InputArray _src, OutputArray _dst,
     _dst.create( map1.size(), src.type() );
     Mat dst = _dst.getMat();
 
-
     CV_Assert( dst.cols < SHRT_MAX && dst.rows < SHRT_MAX && src.cols < SHRT_MAX && src.rows < SHRT_MAX );
 
     if( dst.data == src.data )
@@ -1423,6 +1451,10 @@ void cv::remap( InputArray _src, OutputArray _dst,
     if( interpolation == INTER_AREA )
         interpolation = INTER_LINEAR;
 
+    if (genericWarp(src, Mat(), map1, map2, dst, interpolation, borderType, borderValue, hasRelativeFlag)) {
+        return;
+    }
+
     int type = src.type(), depth = CV_MAT_DEPTH(type);
 
     if (interpolation == INTER_NEAREST && map1.depth() == CV_32F) {
@@ -1575,24 +1607,6 @@ void cv::remap( InputArray _src, OutputArray _dst,
         }
     };
 
-    static RemapFunc cubic_tab[2][CV_DEPTH_MAX] =
-    {
-        {
-            remapBicubic<FixedPtCast<int, uchar, INTER_REMAP_COEF_BITS>, short, INTER_REMAP_COEF_SCALE, false>, 0,
-            remapBicubic<Cast<float, ushort>, float, 1, false>,
-            remapBicubic<Cast<float, short>, float, 1, false>, 0,
-            remapBicubic<Cast<float, float>, float, 1, false>,
-            remapBicubic<Cast<double, double>, float, 1, false>, 0
-        },
-        {
-            remapBicubic<FixedPtCast<int, uchar, INTER_REMAP_COEF_BITS>, short, INTER_REMAP_COEF_SCALE, true>, 0,
-            remapBicubic<Cast<float, ushort>, float, 1, true>,
-            remapBicubic<Cast<float, short>, float, 1, true>, 0,
-            remapBicubic<Cast<float, float>, float, 1, true>,
-            remapBicubic<Cast<double, double>, float, 1, true>, 0
-        }
-    };
-
     static RemapFunc lanczos4_tab[2][8] =
     {
         {
@@ -1621,10 +1635,6 @@ void cv::remap( InputArray _src, OutputArray _dst,
     {
         if( interpolation == INTER_LINEAR )
             ifunc = linear_tab[relativeOptionIndex][depth];
-        else if( interpolation == INTER_CUBIC ){
-            ifunc = cubic_tab[relativeOptionIndex][depth];
-            CV_Assert( _src.channels() <= 4 );
-        }
         else if( interpolation == INTER_LANCZOS4 ){
             ifunc = lanczos4_tab[relativeOptionIndex][depth];
             CV_Assert( _src.channels() <= 4 );
@@ -2052,7 +2062,7 @@ static bool ocl_warpTransform_cols4(InputArray _src, OutputArray _dst, InputArra
 
     if ( !dev.isIntel() || !(type == CV_8UC1) ||
          !(dtype == CV_8UC1) || !(_dst.cols() % 4 == 0) ||
-         (op_type == OCL_OP_PERSPECTIVE && interpolation == INTER_LINEAR && (cn == 1 || cn == 3 || cn == 4)) ||
+         (op_type == OCL_OP_PERSPECTIVE && interpolation != INTER_NEAREST) ||
          !(borderType == cv::BORDER_CONSTANT &&
           (interpolation == cv::INTER_NEAREST || interpolation == cv::INTER_LINEAR || interpolation == cv::INTER_CUBIC)))
         return false;
@@ -2062,7 +2072,7 @@ static bool ocl_warpTransform_cols4(InputArray _src, OutputArray _dst, InputArra
     ocl::ProgramSource program = ocl::imgproc::warp_transform_oclsrc;
     String kernelName = format("warp%s_%s_8u", warp_op[op_type], interpolationMap[interpolation]);
 
-    bool is32f = (interpolation == INTER_CUBIC || interpolation == INTER_LINEAR) && op_type == OCL_OP_AFFINE;
+    bool is32f = interpolation == INTER_CUBIC || interpolation == INTER_LINEAR;
     int wdepth = interpolation == INTER_NEAREST ? depth : std::max(is32f ? CV_32F : CV_32S, depth);
     int sctype = CV_MAKETYPE(wdepth, cn);
 
@@ -2145,9 +2155,7 @@ static bool ocl_warpTransform(InputArray _src, OutputArray _dst, InputArray _M0,
     const char * const kernelName = op_type == OCL_OP_AFFINE ? "warpAffine" : "warpPerspective";
 
     int scalarcn = cn == 3 ? 4 : cn;
-    bool is32f = op_type == OCL_OP_AFFINE ?
-                 /* Affine*/ !dev.isAMD() && (interpolation == INTER_CUBIC || interpolation == INTER_LINEAR) :
-                 /* Perspective*/ interpolation == INTER_LINEAR;
+    bool is32f = interpolation == INTER_CUBIC || interpolation == INTER_LINEAR;
     int wdepth = interpolation == INTER_NEAREST ? depth : std::max(is32f ? CV_32F : CV_32S, depth);
     int sctype = CV_MAKETYPE(wdepth, scalarcn);
 
@@ -2490,6 +2498,10 @@ void cv::warpAffine( InputArray _src, OutputArray _dst,
         M[2] = b1; M[5] = b2;
     }
 
+    if (genericWarp(src, matM, Mat(), Mat(), dst, interpolation, borderType, borderValue, false)) {
+        return;
+    }
+
     hal::warpAffine(src.type(), src.data, src.step, src.cols, src.rows, dst.data, dst.step, dst.cols, dst.rows,
                     M, interpolation, borderType, borderValue.val, hint);
 }
@@ -3009,6 +3021,8 @@ void cv::warpPerspective( InputArray _src, OutputArray _dst, InputArray _M0,
 
     CV_Assert( _src.total() > 0 );
 
+    int interpolation = flags & INTER_MAX;
+
     CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat() &&
                _src.cols() <= SHRT_MAX && _src.rows() <= SHRT_MAX,
                ocl_warpTransform_cols4(_src, _dst, _M0, dsize, flags, borderType, borderValue,
@@ -3027,7 +3041,7 @@ void cv::warpPerspective( InputArray _src, OutputArray _dst, InputArray _M0,
 
     double M[9];
     Mat matM(3, 3, CV_64F, M);
-    int interpolation = flags & INTER_MAX;
+
     if( interpolation == INTER_AREA )
         interpolation = INTER_LINEAR;
 
@@ -3037,6 +3051,10 @@ void cv::warpPerspective( InputArray _src, OutputArray _dst, InputArray _M0,
     if( !(flags & WARP_INVERSE_MAP) )
         invert(matM, matM);
 
+    if (genericWarp(src, matM, Mat(), Mat(), dst, interpolation, borderType, borderValue, false)) {
+        return;
+    }
+
     hal::warpPerspective(src.type(), src.data, src.step, src.cols, src.rows, dst.data, dst.step, dst.cols, dst.rows,
                         matM.ptr<double>(), interpolation, borderType, borderValue.val, hint);
 }

modules/imgproc/src/opencl/remap.cl

@@ -51,6 +51,13 @@
 #endif
 #endif
 
+#define CV_64F 6
+#if defined SRC_DEPTH && SRC_DEPTH == CV_64F
+#define WT1 double
+#else
+#define WT1 float
+#endif
+
 #define noconvert
 
 #if CN != 3
@@ -424,8 +431,8 @@ __kernel void remap_2_32FC1(__global const uchar * srcptr, int src_step, int src
             int sx = convert_int_rtn(X0);
             int sy = convert_int_rtn(Y0);
 
-            float ax = X0 - (float) sx;
-            float ay = Y0 - (float) sy;
+            WT1 ax = (WT1)(X0 - (float) sx);
+            WT1 ay = (WT1)(Y0 - (float) sy);
 
             int2 map_data0 = (int2)(sx, sy);
             int2 map_data1 = (int2)(sx+1, sy);
@@ -493,8 +500,8 @@ __kernel void remap_32FC2(__global const uchar * srcptr, int src_step, int src_o
             int sx = convert_int_rtn(X0);
             int sy = convert_int_rtn(Y0);
 
-            float ax = X0 - (float) sx;
-            float ay = Y0 - (float) sy;
+            WT1 ax = (WT1)(X0 - (float) sx);
+            WT1 ay = (WT1)(Y0 - (float) sy);
 
             int2 map_data0 = (int2)(sx, sy);
             int2 map_data1 = (int2)(sx+1, sy);

modules/imgproc/src/opencl/warp_affine.cl

@@ -51,6 +51,13 @@
 #endif
 #endif
 
+#define CV_64F 6
+#if defined SRC_DEPTH && SRC_DEPTH == CV_64F
+#define WT1 double
+#else
+#define WT1 float
+#endif
+
 #define INTER_BITS 5
 #define INTER_TAB_SIZE (1 << INTER_BITS)
 #define INTER_SCALE 1.f/INTER_TAB_SIZE
@@ -137,8 +144,8 @@ __kernel void warpAffine(__global const uchar * srcptr, int src_step, int src_of
             int sx = convert_short_rtn(X0);
             int sy = convert_short_rtn(Y0);
 
-            float ax = X0 - (CT)sx;
-            float ay = Y0 - (CT)sy;
+            WT1 ax = (WT1)(X0 - (CT)sx);
+            WT1 ay = (WT1)(Y0 - (CT)sy);
 
             WT v0 = scalar, v1 = scalar, v2 = scalar, v3 = scalar;
             if (sx >= 0 && sx < src_cols)
@@ -168,35 +175,16 @@ __kernel void warpAffine(__global const uchar * srcptr, int src_step, int src_of
 
 #elif defined INTER_CUBIC
 
-#ifdef AMD_DEVICE
-
-inline void interpolateCubic( float x, float* coeffs )
+inline void interpolateCubic( float x, WT1* coeffs )
 {
     const float A = -0.75f;
 
-    coeffs[0] = fma(fma(fma(A, (x + 1.f), - 5.0f*A), (x + 1.f), 8.0f*A), x + 1.f, - 4.0f*A);
-    coeffs[1] = fma(fma(A + 2.f, x, - (A + 3.f)), x*x, 1.f);
-    coeffs[2] = fma(fma(A + 2.f, 1.f - x, - (A + 3.f)), (1.f - x)*(1.f - x), 1.f);
-    coeffs[3] = 1.f - coeffs[0] - coeffs[1] - coeffs[2];
+    coeffs[0] = (WT1)fma(fma(fma(A, (x + 1.f), - 5.0f*A), (x + 1.f), 8.0f*A), x + 1.f, - 4.0f*A);
+    coeffs[1] = (WT1)fma(fma(A + 2.f, x, - (A + 3.f)), x*x, 1.f);
+    coeffs[2] = (WT1)fma(fma(A + 2.f, 1.f - x, - (A + 3.f)), (1.f - x)*(1.f - x), 1.f);
+    coeffs[3] = (WT1)(1. - coeffs[0] - coeffs[1] - coeffs[2]);
 }
 
-#else
-
-__constant float coeffs[128] =
-    { 0.000000f, 1.000000f, 0.000000f, 0.000000f, -0.021996f, 0.997841f, 0.024864f, -0.000710f, -0.041199f, 0.991516f, 0.052429f, -0.002747f,
-    -0.057747f, 0.981255f, 0.082466f, -0.005974f, -0.071777f, 0.967285f, 0.114746f, -0.010254f, -0.083427f, 0.949837f, 0.149040f, -0.015450f,
-    -0.092834f, 0.929138f, 0.185120f, -0.021423f, -0.100136f, 0.905418f, 0.222755f, -0.028038f, -0.105469f, 0.878906f, 0.261719f, -0.035156f,
-    -0.108971f, 0.849831f, 0.301781f, -0.042641f, -0.110779f, 0.818420f, 0.342712f, -0.050354f, -0.111031f, 0.784904f, 0.384285f, -0.058159f,
-    -0.109863f, 0.749512f, 0.426270f, -0.065918f, -0.107414f, 0.712471f, 0.468437f, -0.073494f, -0.103821f, 0.674011f, 0.510559f, -0.080750f,
-    -0.099220f, 0.634361f, 0.552406f, -0.087547f, -0.093750f, 0.593750f, 0.593750f, -0.093750f, -0.087547f, 0.552406f, 0.634361f, -0.099220f,
-    -0.080750f, 0.510559f, 0.674011f, -0.103821f, -0.073494f, 0.468437f, 0.712471f, -0.107414f, -0.065918f, 0.426270f, 0.749512f, -0.109863f,
-    -0.058159f, 0.384285f, 0.784904f, -0.111031f, -0.050354f, 0.342712f, 0.818420f, -0.110779f, -0.042641f, 0.301781f, 0.849831f, -0.108971f,
-    -0.035156f, 0.261719f, 0.878906f, -0.105469f, -0.028038f, 0.222755f, 0.905418f, -0.100136f, -0.021423f, 0.185120f, 0.929138f, -0.092834f,
-    -0.015450f, 0.149040f, 0.949837f, -0.083427f, -0.010254f, 0.114746f, 0.967285f, -0.071777f, -0.005974f, 0.082466f, 0.981255f, -0.057747f,
-    -0.002747f, 0.052429f, 0.991516f, -0.041199f, -0.000710f, 0.024864f, 0.997841f, -0.021996f };
-
-#endif
-
 __kernel void warpAffine(__global const uchar * srcptr, int src_step, int src_offset, int src_rows, int src_cols,
                          __global uchar * dstptr, int dst_step, int dst_offset, int dst_rows, int dst_cols,
                          __constant CT * M, ST scalar_)
@@ -206,103 +194,60 @@ __kernel void warpAffine(__global const uchar * srcptr, int src_step, int src_of
 
     if (dx < dst_cols && dy < dst_rows)
     {
-        int tmp = (dx << AB_BITS);
-        int X0 = rint(M[0] * tmp) + rint(fma(M[1], (CT)dy, M[2]) * AB_SCALE) + ROUND_DELTA;
-        int Y0 = rint(M[3] * tmp) + rint(fma(M[4], (CT)dy, M[5]) * AB_SCALE) + ROUND_DELTA;
+        float X0 = (float)fma(M[1], (CT)dy, fma(M[0], (CT)dx, M[2]));
+        float Y0 = (float)fma(M[4], (CT)dy, fma(M[3], (CT)dx, M[5]));
 
-        X0 = X0 >> (AB_BITS - INTER_BITS);
-        Y0 = Y0 >> (AB_BITS - INTER_BITS);
+        int sx = convert_short_rtn(X0);
+        int sy = convert_short_rtn(Y0);
+        float ax = X0 - (float)sx;
+        float ay = Y0 - (float)sy;
 
-        int sx = (short)(X0 >> INTER_BITS) - 1, sy = (short)(Y0 >> INTER_BITS) - 1;
-        int ay = (short)(Y0 & (INTER_TAB_SIZE - 1)), ax = (short)(X0 & (INTER_TAB_SIZE - 1));
+        sx--;
+        sy--;
 
-#ifdef AMD_DEVICE
-        WT v[16];
-        #pragma unroll
-        for (int y = 0; y < 4; y++)
-        {
-            if (sy+y >= 0 && sy+y < src_rows)
-            {
-                #pragma unroll
-                for (int x = 0; x < 4; x++)
-                    v[mad24(y, 4, x)] = sx+x >= 0 && sx+x < src_cols ?
-                        CONVERT_TO_WT(loadpix(srcptr + mad24(sy+y, src_step, mad24(sx+x, pixsize, src_offset)))) : scalar;
-            }
-            else
-            {
-                #pragma unroll
-                for (int x = 0; x < 4; x++)
-                    v[mad24(y, 4, x)] = scalar;
-            }
-        }
-
-        float tab1y[4], tab1x[4];
-
-        float ayy = INTER_SCALE * ay;
-        float axx = INTER_SCALE * ax;
-        interpolateCubic(ayy, tab1y);
-        interpolateCubic(axx, tab1x);
-
-        int dst_index = mad24(dy, dst_step, mad24(dx, pixsize, dst_offset));
+        WT1 taby[4], tabx[4];
+        interpolateCubic(ay, taby);
+        interpolateCubic(ax, tabx);
 
         WT sum = (WT)(0);
-#if SRC_DEPTH <= 4
-        int itab[16];
 
-        #pragma unroll
-        for (int i = 0; i < 16; i++)
-            itab[i] = rint(tab1y[(i>>2)] * tab1x[(i&3)] * INTER_REMAP_COEF_SCALE);
-
-        #pragma unroll
-        for (int i = 0; i < 16; i++)
-            sum = mad24(v[i], itab[i], sum);
-        storepix(CONVERT_TO_T( (sum + (1 << (INTER_REMAP_COEF_BITS-1))) >> INTER_REMAP_COEF_BITS ), dstptr + dst_index);
-#else
-        #pragma unroll
-        for (int i = 0; i < 16; i++)
-            sum = fma(v[i], tab1y[(i>>2)] * tab1x[(i&3)], sum);
-        storepix(CONVERT_TO_T( sum ), dstptr + dst_index);
-#endif
-#else // INTEL_DEVICE
-        __constant float * coeffs_y = coeffs + (ay << 2), * coeffs_x = coeffs + (ax << 2);
-
-        int src_index0 = mad24(sy, src_step, mad24(sx, pixsize, src_offset)), src_index;
-        int dst_index = mad24(dy, dst_step, mad24(dx, pixsize, dst_offset));
-
-        WT sum = (WT)(0), xsum;
-        #pragma unroll
-        for (int y = 0; y < 4; y++)
-        {
-            src_index = mad24(y, src_step, src_index0);
-            if (sy + y >= 0 && sy + y < src_rows)
-            {
-                xsum = (WT)(0);
-                if (sx >= 0 && sx + 4 < src_cols)
-                {
-#if SRC_DEPTH == 0 && CN == 1
-                    uchar4 value = vload4(0, srcptr + src_index);
-                    xsum = dot(convert_float4(value), (float4)(coeffs_x[0], coeffs_x[1], coeffs_x[2], coeffs_x[3]));
-#else
-                    #pragma unroll
-                    for (int x = 0; x < 4; x++)
-                        xsum = fma(CONVERT_TO_WT(loadpix(srcptr + mad24(x, pixsize, src_index))), coeffs_x[x], xsum);
-#endif
-                }
-                else
-                {
-                    #pragma unroll
-                    for (int x = 0; x < 4; x++)
-                        xsum = fma(sx + x >= 0 && sx + x < src_cols ?
-                                   CONVERT_TO_WT(loadpix(srcptr + mad24(x, pixsize, src_index))) : scalar, coeffs_x[x], xsum);
-                }
-                sum = fma(xsum, coeffs_y[y], sum);
+        if (0 <= sy && sy + 4 <= src_rows &&
+            0 <= sx && sx + 4 <= src_cols) {
+            #pragma unroll
+            for (int y = 0; y < 4; y++) {
+                int row_offset = mad24(sy+y, src_step, src_offset);
+                WT v0 = CONVERT_TO_WT(loadpix(srcptr + mad24(sx, pixsize, row_offset)));
+                WT v1 = CONVERT_TO_WT(loadpix(srcptr + mad24(sx + 1, pixsize, row_offset)));
+                WT v2 = CONVERT_TO_WT(loadpix(srcptr + mad24(sx + 2, pixsize, row_offset)));
+                WT v3 = CONVERT_TO_WT(loadpix(srcptr + mad24(sx + 3, pixsize, row_offset)));
+                WT wsum = (WT)(0);
+                wsum = fma(v0, tabx[0], wsum);
+                wsum = fma(v1, tabx[1], wsum);
+                wsum = fma(v2, tabx[2], wsum);
+                wsum = fma(v3, tabx[3], wsum);
+                sum = fma(wsum, taby[y], sum);
+            }
+        }
+        else {
+            #pragma unroll
+            for (int y = 0; y < 4; y++) {
+                if (sy+y >= 0 && sy+y < src_rows) {
+                    int row_offset = mad24(sy+y, src_step, src_offset);
+                    #pragma unroll
+                    for (int x = 0; x < 4; x++) {
+                        WT v = sx+x >= 0 && sx+x < src_cols ?
+                            CONVERT_TO_WT(loadpix(srcptr + mad24(sx + x, pixsize, row_offset))) : scalar;
+                        sum = fma(v, taby[y] * tabx[x], sum);
+                    }
+                }
+                else {
+                    sum = fma(scalar, taby[y], sum);
+                }
             }
-            else
-                sum = fma(scalar, coeffs_y[y], sum);
         }
 
-        storepix(CONVERT_TO_T(sum), dstptr + dst_index);
-#endif
+        int dst_index = mad24(dy, dst_step, mad24(dx, pixsize, dst_offset));
+        storepix(CONVERT_TO_T( sum ), dstptr + dst_index);
     }
 }
 

modules/imgproc/src/opencl/warp_perspective.cl

@@ -51,6 +51,13 @@
 #endif
 #endif
 
+#define CV_64F 6
+#if defined SRC_DEPTH && SRC_DEPTH == CV_64F
+#define WT1 double
+#else
+#define WT1 float
+#endif
+
 #define INTER_BITS 5
 #define INTER_TAB_SIZE (1 << INTER_BITS)
 #define INTER_SCALE 1.f / INTER_TAB_SIZE
@@ -150,14 +157,14 @@ __kernel void warpPerspective(__global const uchar * srcptr, int src_step, int s
 
 #elif defined INTER_CUBIC
 
-inline void interpolateCubic( float x, float* coeffs )
+inline void interpolateCubic( float x, WT1* coeffs )
 {
     const float A = -0.75f;
 
-    coeffs[0] = ((A*(x + 1.f) - 5.0f*A)*(x + 1.f) + 8.0f*A)*(x + 1.f) - 4.0f*A;
-    coeffs[1] = ((A + 2.f)*x - (A + 3.f))*x*x + 1.f;
-    coeffs[2] = ((A + 2.f)*(1.f - x) - (A + 3.f))*(1.f - x)*(1.f - x) + 1.f;
-    coeffs[3] = 1.f - coeffs[0] - coeffs[1] - coeffs[2];
+    coeffs[0] = (WT1)fma(fma(fma(A, (x + 1.f), - 5.0f*A), (x + 1.f), 8.0f*A), x + 1.f, - 4.0f*A);
+    coeffs[1] = (WT1)fma(fma(A + 2.f, x, - (A + 3.f)), x*x, 1.f);
+    coeffs[2] = (WT1)fma(fma(A + 2.f, 1.f - x, - (A + 3.f)), (1.f - x)*(1.f - x), 1.f);
+    coeffs[3] = (WT1)(1. - coeffs[0] - coeffs[1] - coeffs[2]);
 }
 
 __kernel void warpPerspective(__global const uchar * srcptr, int src_step, int src_offset, int src_rows, int src_cols,
@@ -169,52 +176,61 @@ __kernel void warpPerspective(__global const uchar * srcptr, int src_step, int s
 
     if (dx < dst_cols && dy < dst_rows)
     {
-        CT X0 = M[0] * dx + M[1] * dy + M[2];
-        CT Y0 = M[3] * dx + M[4] * dy + M[5];
-        CT W = M[6] * dx + M[7] * dy + M[8];
-        W = W != 0.0f ? INTER_TAB_SIZE / W : 0.0f;
-        int X = rint(X0 * W), Y = rint(Y0 * W);
+        float W = (float)fma(M[7], (CT)dy, fma(M[6], (CT)dx, M[8]));
+        float X0 = (float)fma(M[1], (CT)dy, fma(M[0], (CT)dx, M[2])) / W;
+        float Y0 = (float)fma(M[4], (CT)dy, fma(M[3], (CT)dx, M[5])) / W;
 
-        short sx = convert_short_sat(X >> INTER_BITS) - 1;
-        short sy = convert_short_sat(Y >> INTER_BITS) - 1;
-        short ay = (short)(Y & (INTER_TAB_SIZE-1));
-        short ax = (short)(X & (INTER_TAB_SIZE-1));
+        int sx = convert_short_rtn(X0);
+        int sy = convert_short_rtn(Y0);
+        float ax = X0 - (float)sx;
+        float ay = Y0 - (float)sy;
 
-        WT v[16];
-        #pragma unroll
-        for (int y = 0; y < 4; y++)
-            #pragma unroll
-            for (int x = 0; x < 4; x++)
-                v[mad24(y, 4, x)] = (sx+x >= 0 && sx+x < src_cols && sy+y >= 0 && sy+y < src_rows) ?
-                    CONVERT_TO_WT(loadpix(srcptr + mad24(sy+y, src_step, src_offset + (sx+x) * pixsize))) : scalar;
+        sx--;
+        sy--;
 
-        float tab1y[4], tab1x[4];
-
-        float ayy = INTER_SCALE * ay;
-        float axx = INTER_SCALE * ax;
-        interpolateCubic(ayy, tab1y);
-        interpolateCubic(axx, tab1x);
-
-        int dst_index = mad24(dy, dst_step, dst_offset + dx * pixsize);
+        WT1 taby[4], tabx[4];
+        interpolateCubic(ay, taby);
+        interpolateCubic(ax, tabx);
 
         WT sum = (WT)(0);
-#if SRC_DEPTH <= 4
-        int itab[16];
 
-        #pragma unroll
-        for (int i = 0; i < 16; i++)
-            itab[i] = rint(tab1y[(i>>2)] * tab1x[(i&3)] * INTER_REMAP_COEF_SCALE);
+        if (0 <= sy && sy + 4 <= src_rows &&
+            0 <= sx && sx + 4 <= src_cols) {
+            #pragma unroll
+            for (int y = 0; y < 4; y++) {
+                int row_offset = mad24(sy+y, src_step, src_offset);
+                WT v0 = CONVERT_TO_WT(loadpix(srcptr + mad24(sx, pixsize, row_offset)));
+                WT v1 = CONVERT_TO_WT(loadpix(srcptr + mad24(sx + 1, pixsize, row_offset)));
+                WT v2 = CONVERT_TO_WT(loadpix(srcptr + mad24(sx + 2, pixsize, row_offset)));
+                WT v3 = CONVERT_TO_WT(loadpix(srcptr + mad24(sx + 3, pixsize, row_offset)));
+                WT wsum = (WT)(0);
+                wsum = fma(v0, tabx[0], wsum);
+                wsum = fma(v1, tabx[1], wsum);
+                wsum = fma(v2, tabx[2], wsum);
+                wsum = fma(v3, tabx[3], wsum);
+                sum = fma(wsum, taby[y], sum);
+            }
+        }
+        else {
+            #pragma unroll
+            for (int y = 0; y < 4; y++) {
+                if (sy+y >= 0 && sy+y < src_rows) {
+                    int row_offset = mad24(sy+y, src_step, src_offset);
+                    #pragma unroll
+                    for (int x = 0; x < 4; x++) {
+                        WT v = sx+x >= 0 && sx+x < src_cols ?
+                            CONVERT_TO_WT(loadpix(srcptr + mad24(sx + x, pixsize, row_offset))) : scalar;
+                        sum = fma(v, taby[y] * tabx[x], sum);
+                    }
+                }
+                else {
+                    sum = fma(scalar, taby[y], sum);
+                }
+            }
+        }
 
-        #pragma unroll
-        for (int i = 0; i < 16; i++)
-            sum += v[i] * itab[i];
-        storepix(CONVERT_TO_T( (sum + (1 << (INTER_REMAP_COEF_BITS-1))) >> INTER_REMAP_COEF_BITS ), dstptr + dst_index);
-#else
-        #pragma unroll
-        for (int i = 0; i < 16; i++)
-            sum += v[i] * tab1y[(i>>2)] * tab1x[(i&3)];
+        int dst_index = mad24(dy, dst_step, mad24(dx, pixsize, dst_offset));
         storepix(CONVERT_TO_T( sum ), dstptr + dst_index);
-#endif
     }
 }
 

modules/imgproc/src/warp_common.hpp

@@ -8,4 +8,15 @@
 #include "warp_common.vector.hpp"
 #include "warp_common.scalar.hpp"
 
+namespace cv {
+
+typedef void (*ImgWarpFunc)(const float* x, const float* y, int len,
+                            const void* src, size_t srcstep, Size size,
+                            void* dst, const float* params,
+                            int borderType, const void* borderVal);
+
+ImgWarpFunc getImgWarpFunc(int type, int interpolation);
+
+}
+
 #endif // __OPENCV_IMGPROC_WARP_COMMON_HPP__

modules/imgproc/test/ocl/test_warp.cpp

@@ -107,7 +107,7 @@ PARAM_TEST_CASE(WarpTestBase, MatType, Interpolation, bool, bool)
     void Near(double threshold = 0.0)
     {
         if (depth < CV_32F)
-            EXPECT_MAT_N_DIFF(dst_roi, udst_roi, cvRound(dst_roi.total()*threshold));
+            EXPECT_MAT_N_DIFF_EPS(dst_roi, udst_roi, 1, cvRound(dst_roi.total()*threshold));
         else
             OCL_EXPECT_MATS_NEAR_RELATIVE(dst, threshold);
     }

modules/imgproc/test/test_imgwarp.cpp

@@ -1114,5 +1114,90 @@ TEST(Imgproc_getPerspectiveTransform, issue_26916)
     EXPECT_MAT_NEAR(obtained_homogeneous_dst_points, expected_homogeneous_dst_points, 1e-10);
 }
 
+static void rotation2affine(float scale, float angle, float cx, float cy, float* M)
+{
+    // (x - cx)*cos(a) + (y - cy)*sin(a) + cx
+    // -(x - cx)*sin(a) + (y - cy)*cos(a) + cy
+    float ca = cosf(angle), sa = sinf(angle);
+    M[0] = scale*ca; M[1] = scale*sa; M[2] = scale*(-cx*ca - cy*sa) + cx;
+    M[3] = -scale*sa; M[4] = scale*ca; M[5] = scale*(cx*sa - cy*ca) + cy;
+}
+
+TEST(Imgproc_Warping, DISABLED_playground)
+{
+    int imgtype = CV_32F;
+    int imgcn = 3;
+    bool useOpenCL = true;
+
+    auto ts = cvtest::TS::ptr();
+    Mat img0 = imread(string(ts->get_data_path()) + "stereomatching/datasets/tsukuba/im2.png", 1), img1, img;
+    int iangle = -1;
+    int borderType = BORDER_CONSTANT;
+    Scalar borderValue(0, 128, 0);
+
+    double cvtscale = imgtype == CV_16U ? 256. : imgtype == CV_32F ? 1./255 : 1.;
+    if (imgcn == 1) {
+        cvtColor(img0, img1, COLOR_BGR2GRAY);
+    } else if (imgcn == 4) {
+        cvtColor(img0, img1, COLOR_BGR2BGRA);
+    } else if (imgcn == 3) {
+        img1 = img0;
+    } else {
+        CV_Assert(imgcn == 2);
+        std::vector<Mat> ch;
+        split(img0, ch);
+        ch.pop_back();
+        merge(ch, img1);
+    }
+    img1.convertTo(img, imgtype, cvtscale);
+    Mat canvas0(img.size(), imgtype), canvas8;
+    float cx = img.cols*0.5f, cy = img.rows*0.5f;
+    if (img.depth() == CV_32F) {
+        borderValue = Scalar(100*cvtscale, 0*cvtscale, 100*cvtscale);
+    }
+    UMat uimg, ucanvas;
+    if (useOpenCL) {
+        img.copyTo(uimg);
+    }
+
+    for(;;) {
+        iangle = (iangle + 1) % (360*4);
+        float angle = float(iangle*CV_PI/180.f*0.25f);
+        float scale = float(1 + 0.2f*sin(angle));
+        float Mdata[9] = {1, 0, 0, 0, 1, 0, 0, 0, 1};
+        rotation2affine(scale, angle, cx, cy, Mdata);
+        Mat M(2, 3, CV_32F, Mdata);
+
+        double t0 = (double)getTickCount();
+        if (!useOpenCL) {
+            warpAffine(img, canvas0, M, canvas0.size(), INTER_CUBIC, borderType, borderValue);
+        } else {
+            warpAffine(uimg, ucanvas, M, uimg.size(), INTER_CUBIC, borderType, borderValue);
+            ocl::finish();
+        }
+        t0 = (double)getTickCount() - t0;
+
+        if (useOpenCL) {
+            ucanvas.copyTo(canvas0);
+        }
+        canvas0.convertTo(canvas8, CV_8U, 1./cvtscale);
+        if (canvas8.channels() == 2) {
+            std::vector<Mat> ch;
+            split(canvas8, ch);
+            ch.push_back(Mat::zeros(canvas8.size(), CV_8U));
+            merge(ch, canvas8);
+        }
+        printf("opencv time = %.1fms\n", t0*1000./getTickFrequency());
+        imshow("result (opencv)", canvas8);
+        int c = waitKey(1);
+        if (c < 0)
+            continue;
+        if ((c & 255) == 27)
+            break;
+        if ((waitKey() & 255) == 27)
+            break;
+    }
+}
+
 }} // namespace
 /* End of file. */

modules/imgproc/test/test_imgwarp_strict.cpp

@@ -231,13 +231,13 @@ void CV_ImageWarpBaseTest::run(int)
 float CV_ImageWarpBaseTest::get_success_error_level(int _interpolation, int) const
 {
     if (_interpolation == INTER_CUBIC)
-        return 1.0f;
+        return 2.0f;
     else if (_interpolation == INTER_LANCZOS4)
         return 1.0f;
     else if (_interpolation == INTER_NEAREST)
         return 255.0f;  // FIXIT: check is not reliable for Black/White (0/255) images
     else if (_interpolation == INTER_AREA)
-        return 2.0f;
+        return 1.0f;
     else
         return 1.0f;
 }
@@ -304,7 +304,7 @@ void CV_ImageWarpBaseTest::validate_results() const
 
                 const int radius = 3;
                 int rmin = MAX(dy - radius, 0), rmax = MIN(dy + radius, dsize.height);
-                int cmin = MAX(dx / cn - radius, 0), cmax = MIN(dx / cn + radius, dsize.width);
+                int cmin = MAX(dx / cn - radius, 0), cmax = MIN(dx / cn + radius, dsize.width/cn);
 
                 std::cout << "opencv result:\n" << dst(Range(rmin, rmax), Range(cmin, cmax)) << std::endl;
                 std::cout << "reference result:\n" << reference_dst(Range(rmin, rmax), Range(cmin, cmax)) << std::endl;
@@ -851,16 +851,23 @@ void CV_Remap_Test::run_func()
 
 void CV_Remap_Test::convert_maps()
 {
-    if (mapx.type() != CV_16SC2)
-        convertMaps(mapx.clone(), mapy.clone(), mapx, mapy, CV_16SC2, interpolation == INTER_NEAREST);
-    else if (interpolation != INTER_NEAREST)
-        if (mapy.type() != CV_16UC1)
-            mapy.clone().convertTo(mapy, CV_16UC1);
+    if (interpolation == INTER_CUBIC) {
+        if (mapx.type() != CV_32FC1)
+            convertMaps(mapx.clone(), mapy.clone(), mapx, mapy, CV_32FC1, interpolation == INTER_NEAREST);
+    }
+    else {
+        if (mapx.type() != CV_16SC2)
+            convertMaps(mapx.clone(), mapy.clone(), mapx, mapy, CV_16SC2, interpolation == INTER_NEAREST);
+        else if (interpolation != INTER_NEAREST) {
+            if (mapy.type() != CV_16UC1)
+                mapy.clone().convertTo(mapy, CV_16UC1);
+        }
 
-    if (interpolation == INTER_NEAREST)
-        mapy = Mat();
-    CV_Assert(((interpolation == INTER_NEAREST && mapy.empty()) || mapy.type() == CV_16UC1 ||
-               mapy.type() == CV_16SC1) && mapx.type() == CV_16SC2);
+        if (interpolation == INTER_NEAREST)
+            mapy = Mat();
+        CV_Assert(((interpolation == INTER_NEAREST && mapy.empty()) || mapy.type() == CV_16UC1 ||
+                   mapy.type() == CV_16SC1) && mapx.type() == CV_16SC2);
+    }
 }
 
 const char* CV_Remap_Test::borderType_to_string() const
@@ -1164,8 +1171,10 @@ void CV_Remap_Test::remap_nearest(const Mat& _src, Mat& _dst)
 
 void CV_Remap_Test::remap_generic(const Mat& _src, Mat& _dst)
 {
-    CV_Assert(mapx.type() == CV_16SC2 && mapy.type() == CV_16UC1);
+    CV_Assert((mapx.type() == CV_16SC2 && mapy.type() == CV_16UC1) ||
+              (mapx.type() == CV_32FC1 && mapy.type() == CV_32FC1));
 
+    bool fixedpt = mapx.type() == CV_16SC2;
     int ksize = 2;
     if (interpolation == INTER_CUBIC)
         ksize = 4;
@@ -1185,19 +1194,42 @@ void CV_Remap_Test::remap_generic(const Mat& _src, Mat& _dst)
 
     for (int dy = 0; dy < dsize.height; ++dy)
     {
-        const short* yMx = mapx.ptr<short>(dy);
-        const ushort* yMy = mapy.ptr<ushort>(dy);
+        const short* Mxy = nullptr;
+        const ushort* Mab = nullptr;
+        const float* Mx = nullptr;
+        const float* My = nullptr;
+
+        if (fixedpt) {
+            Mxy = mapx.ptr<short>(dy);
+            Mab = mapy.ptr<ushort>(dy);
+        } else {
+            Mx = mapx.ptr<float>(dy);
+            My = mapy.ptr<float>(dy);
+        }
 
         float* yD = _dst.ptr<float>(dy);
 
         for (int dx = 0; dx < dsize.width; ++dx)
         {
             float* xyD = yD + dx * cn;
-            float sx = yMx[dx * 2], sy = yMx[dx * 2 + 1];
-            int isx = cvFloor(sx), isy = cvFloor(sy);
+            int isx, isy;
+            float alpha, beta;
+            if (fixedpt) {
+                isx = Mxy[dx * 2];
+                isy = Mxy[dx * 2 + 1];
+                alpha = (Mab[dx] & (INTER_TAB_SIZE - 1)) / static_cast<float>(INTER_TAB_SIZE);
+                beta = ((Mab[dx] >> INTER_BITS) & (INTER_TAB_SIZE - 1)) / static_cast<float>(INTER_TAB_SIZE);
+            } else {
+                float sx = Mx[dx];
+                float sy = My[dx];
+                isx = cvFloor(sx);
+                isy = cvFloor(sy);
+                alpha = sx - isx;
+                beta = sy - isy;
+            }
 
-            inter_func((yMy[dx] & (INTER_TAB_SIZE - 1)) / static_cast<float>(INTER_TAB_SIZE), w);
-            inter_func(((yMy[dx] >> INTER_BITS) & (INTER_TAB_SIZE - 1)) / static_cast<float>(INTER_TAB_SIZE), w + ksize);
+            inter_func(alpha, w);
+            inter_func(beta, w + ksize);
 
             isx -= ofs;
             isy -= ofs;
@@ -1332,7 +1364,7 @@ void CV_WarpAffine_Test::run_func()
 
 float CV_WarpAffine_Test::get_success_error_level(int _interpolation, int _depth) const
 {
-    return _depth == CV_8U ? 0.f : CV_ImageWarpBaseTest::get_success_error_level(_interpolation, _depth);
+    return _depth == CV_8U ? 2.f : CV_ImageWarpBaseTest::get_success_error_level(_interpolation, _depth);
 }
 
 void CV_WarpAffine_Test::run_reference_func()
@@ -1390,6 +1422,44 @@ void CV_WarpAffine_Test::new_warpAffine(const Mat &_src, Mat &_dst, const Mat &t
     }
 }
 
+static void computeWarpMaps(const Mat& M, Mat& mapx, Mat& mapy)
+{
+    CV_Assert(M.type() == CV_64F);
+    CV_Assert(M.size() == Size(3, 2) ||
+              M.size() == Size(3, 3));
+    CV_Assert(mapx.size() == mapy.size());
+    CV_Assert(mapx.type() == CV_32FC1 && mapy.type() == CV_32FC1);
+
+    bool perspective = M.rows == 3;
+    const double* dataM = M.ptr<double>();
+    Size dsize = mapx.size();
+
+    for (int dy = 0; dy < dsize.height; ++dy)
+    {
+        float* mapxrow = mapx.ptr<float>(dy);
+        float* mapyrow = mapy.ptr<float>(dy);
+        for (int dx = 0; dx < dsize.width; ++dx)
+        {
+            double x = dataM[0] * dx + dataM[1] * dy + dataM[2];
+            double y = dataM[3] * dx + dataM[4] * dy + dataM[5];
+
+            if (perspective) {
+                double z = dataM[6] * dx + dataM[7] * dy + dataM[8];
+                if (z != 0.0) {
+                    x /= z;
+                    y /= z;
+                } else {
+                    x = 0.0;
+                    y = 0.0;
+                }
+            }
+
+            mapxrow[dx] = float(x);
+            mapyrow[dx] = float(y);
+        }
+    }
+}
+
 void CV_WarpAffine_Test::warpAffine(const Mat& _src, Mat& _dst)
 {
     Size dsize = _dst.size();
@@ -1405,12 +1475,6 @@ void CV_WarpAffine_Test::warpAffine(const Mat& _src, Mat& _dst)
     if (inter == INTER_AREA)
         inter = INTER_LINEAR;
 
-    mapx.create(dsize, CV_16SC2);
-    if (inter != INTER_NEAREST)
-        mapy.create(dsize, CV_16SC1);
-    else
-        mapy = Mat();
-
     if (!(interpolation & cv::WARP_INVERSE_MAP))
         invertAffineTransform(tM.clone(), tM);
 
@@ -1425,32 +1489,44 @@ void CV_WarpAffine_Test::warpAffine(const Mat& _src, Mat& _dst)
         }
     }
 
-    const int AB_BITS = MAX(10, (int)INTER_BITS);
-    const int AB_SCALE = 1 << AB_BITS;
-    int round_delta = (inter == INTER_NEAREST) ? AB_SCALE / 2 : (AB_SCALE / INTER_TAB_SIZE / 2);
+    if (inter == INTER_CUBIC) {
+        mapx.create(dsize, CV_32FC1);
+        mapy.create(dsize, CV_32FC1);
+        computeWarpMaps(tM, mapx, mapy);
+    } else {
+        mapx.create(dsize, CV_16SC2);
+        if (inter != INTER_NEAREST)
+            mapy.create(dsize, CV_16SC1);
+        else
+            mapy = Mat();
 
-    const softdouble* data_tM = tM.ptr<softdouble>(0);
-    for (int dy = 0; dy < dsize.height; ++dy)
-    {
-        short* yM = mapx.ptr<short>(dy);
-        for (int dx = 0; dx < dsize.width; ++dx, yM += 2)
+        const int AB_BITS = MAX(10, (int)INTER_BITS);
+        const int AB_SCALE = 1 << AB_BITS;
+        int round_delta = (inter == INTER_NEAREST) ? AB_SCALE / 2 : (AB_SCALE / INTER_TAB_SIZE / 2);
+
+        const softdouble* data_tM = tM.ptr<softdouble>(0);
+        for (int dy = 0; dy < dsize.height; ++dy)
         {
-            int v1 = saturate_cast<int>(saturate_cast<int>(data_tM[0] * dx * AB_SCALE) +
-                    saturate_cast<int>((data_tM[1] * dy + data_tM[2]) * AB_SCALE) + round_delta),
+            short* yM = mapx.ptr<short>(dy);
+            for (int dx = 0; dx < dsize.width; ++dx, yM += 2)
+            {
+                int v1 = saturate_cast<int>(saturate_cast<int>(data_tM[0] * dx * AB_SCALE) +
+                                            saturate_cast<int>((data_tM[1] * dy + data_tM[2]) * AB_SCALE) + round_delta),
                 v2 = saturate_cast<int>(saturate_cast<int>(data_tM[3] * dx * AB_SCALE) +
-                    saturate_cast<int>((data_tM[4] * dy + data_tM[5]) * AB_SCALE) + round_delta);
-            v1 >>= AB_BITS - INTER_BITS;
-            v2 >>= AB_BITS - INTER_BITS;
+                                        saturate_cast<int>((data_tM[4] * dy + data_tM[5]) * AB_SCALE) + round_delta);
+                v1 >>= AB_BITS - INTER_BITS;
+                v2 >>= AB_BITS - INTER_BITS;
 
-            yM[0] = saturate_cast<short>(v1 >> INTER_BITS);
-            yM[1] = saturate_cast<short>(v2 >> INTER_BITS);
+                yM[0] = saturate_cast<short>(v1 >> INTER_BITS);
+                yM[1] = saturate_cast<short>(v2 >> INTER_BITS);
 
-            if (inter != INTER_NEAREST)
-                mapy.ptr<short>(dy)[dx] = ((v2 & (INTER_TAB_SIZE - 1)) * INTER_TAB_SIZE + (v1 & (INTER_TAB_SIZE - 1)));
+                if (inter != INTER_NEAREST)
+                    mapy.ptr<short>(dy)[dx] = ((v2 & (INTER_TAB_SIZE - 1)) * INTER_TAB_SIZE + (v1 & (INTER_TAB_SIZE - 1)));
+            }
         }
-    }
 
-    CV_Assert(mapx.type() == CV_16SC2 && ((inter == INTER_NEAREST && mapy.empty()) || mapy.type() == CV_16SC1));
+        CV_Assert(mapx.type() == CV_16SC2 && ((inter == INTER_NEAREST && mapy.empty()) || mapy.type() == CV_16SC1));
+    }
     cv::remap(_src, _dst, mapx, mapy, inter, borderType, borderValue);
 }
 
@@ -1612,41 +1688,47 @@ void CV_WarpPerspective_Test::warpPerspective(const Mat& _src, Mat& _dst)
         }
     }
 
-    mapx.create(dsize, CV_16SC2);
-    if (inter != INTER_NEAREST)
-        mapy.create(dsize, CV_16SC1);
-    else
-        mapy = Mat();
+    if (inter == INTER_CUBIC) {
+        mapx.create(dsize, CV_32FC1);
+        mapy.create(dsize, CV_32FC1);
+        computeWarpMaps(M, mapx, mapy);
+    } else {
+        mapx.create(dsize, CV_16SC2);
+        if (inter != INTER_NEAREST)
+            mapy.create(dsize, CV_16SC1);
+        else
+            mapy = Mat();
 
-    double* tM = M.ptr<double>(0);
-    for (int dy = 0; dy < dsize.height; ++dy)
-    {
-        short* yMx = mapx.ptr<short>(dy);
-
-        for (int dx = 0; dx < dsize.width; ++dx, yMx += 2)
+        double* tM = M.ptr<double>(0);
+        for (int dy = 0; dy < dsize.height; ++dy)
         {
-            double den = tM[6] * dx + tM[7] * dy + tM[8];
-            den = den ? 1.0 / den : 0.0;
+            short* yMx = mapx.ptr<short>(dy);
 
-            if (inter == INTER_NEAREST)
+            for (int dx = 0; dx < dsize.width; ++dx, yMx += 2)
             {
-                yMx[0] = saturate_cast<short>((tM[0] * dx + tM[1] * dy + tM[2]) * den);
-                yMx[1] = saturate_cast<short>((tM[3] * dx + tM[4] * dy + tM[5]) * den);
-                continue;
+                double den = tM[6] * dx + tM[7] * dy + tM[8];
+                den = den ? 1.0 / den : 0.0;
+
+                if (inter == INTER_NEAREST)
+                {
+                    yMx[0] = saturate_cast<short>((tM[0] * dx + tM[1] * dy + tM[2]) * den);
+                    yMx[1] = saturate_cast<short>((tM[3] * dx + tM[4] * dy + tM[5]) * den);
+                    continue;
+                }
+
+                den *= static_cast<double>(INTER_TAB_SIZE);
+                int v0 = saturate_cast<int>((tM[0] * dx + tM[1] * dy + tM[2]) * den);
+                int v1 = saturate_cast<int>((tM[3] * dx + tM[4] * dy + tM[5]) * den);
+
+                yMx[0] = saturate_cast<short>(v0 >> INTER_BITS);
+                yMx[1] = saturate_cast<short>(v1 >> INTER_BITS);
+                mapy.ptr<short>(dy)[dx] = saturate_cast<short>((v1 & (INTER_TAB_SIZE - 1)) *
+                                                               INTER_TAB_SIZE + (v0 & (INTER_TAB_SIZE - 1)));
             }
-
-            den *= static_cast<double>(INTER_TAB_SIZE);
-            int v0 = saturate_cast<int>((tM[0] * dx + tM[1] * dy + tM[2]) * den);
-            int v1 = saturate_cast<int>((tM[3] * dx + tM[4] * dy + tM[5]) * den);
-
-            yMx[0] = saturate_cast<short>(v0 >> INTER_BITS);
-            yMx[1] = saturate_cast<short>(v1 >> INTER_BITS);
-            mapy.ptr<short>(dy)[dx] = saturate_cast<short>((v1 & (INTER_TAB_SIZE - 1)) *
-                    INTER_TAB_SIZE + (v0 & (INTER_TAB_SIZE - 1)));
         }
-    }
 
-    CV_Assert(mapx.type() == CV_16SC2 && ((inter == INTER_NEAREST && mapy.empty()) || mapy.type() == CV_16SC1));
+        CV_Assert(mapx.type() == CV_16SC2 && ((inter == INTER_NEAREST && mapy.empty()) || mapy.type() == CV_16SC1));
+    }
     cv::remap(_src, _dst, mapx, mapy, inter, borderType, borderValue);
 }
 

modules/ts/include/opencv2/ts/ocl_test.hpp

@@ -110,6 +110,18 @@ do \
         << "Size: " << (mat1).size() << std::endl; \
 } while ((void)0, 0)
 
+#define EXPECT_MAT_N_DIFF_EPS(mat1, mat2, eps, num) \
+do \
+{ \
+    ASSERT_EQ(mat1.type(), mat2.type()); \
+    ASSERT_EQ(mat1.size(), mat2.size()); \
+    Mat diff, mask; \
+    absdiff(mat1, mat2, diff); \
+    cv::compare(diff.reshape(1), Scalar::all(eps), mask, CMP_GT); \
+    EXPECT_LE(countNonZero(mask), num) \
+    << "Size: " << mat1.size() << std::endl; \
+} while ((void)0, 0)
+
 #define EXPECT_MAT_N_DIFF(mat1, mat2, num) \
 do \
 { \