diff --git a/modules/imgproc/src/opencl/remap.cl b/modules/imgproc/src/opencl/remap.cl index 668ea9519e..1d02e583e3 100644 --- a/modules/imgproc/src/opencl/remap.cl +++ b/modules/imgproc/src/opencl/remap.cl @@ -437,8 +437,9 @@ __kernel void remap_2_32FC1(__global const uchar * srcptr, int src_step, int src #if defined BORDER_CONSTANT float xf = map1[0], yf = map2[0]; - int sx = (convert_int_sat_rtz(mad(xf, (float)INTER_TAB_SIZE, 0.5f)) >> INTER_BITS); - int sy = (convert_int_sat_rtz(mad(yf, (float)INTER_TAB_SIZE, 0.5f)) >> INTER_BITS); + // Coordinates can be negative, so truncation toward zero gives incorrect indices. + int sx = (convert_int_sat_rte(xf * (float)INTER_TAB_SIZE) >> INTER_BITS); + int sy = (convert_int_sat_rte(yf * (float)INTER_TAB_SIZE) >> INTER_BITS); #if WARP_RELATIVE sx += x; sy += y; @@ -464,7 +465,7 @@ __kernel void remap_2_32FC1(__global const uchar * srcptr, int src_step, int src #else #pragma unroll for (int xp = 0; xp < 2; ++xp) - xsum = fma(CONVERT_TO_WT(loadpix(srcptr + mad24(xp, TSIZE, src_index))), coeffs_x[xp], xsum); + xsum = fma(CONVERT_TO_WT(loadpix(srcptr + mad24(xp, TSIZE, src_index))), (WT)(coeffs_x[xp]), xsum); #endif } else @@ -472,12 +473,12 @@ __kernel void remap_2_32FC1(__global const uchar * srcptr, int src_step, int src #pragma unroll for (int xp = 0; xp < 2; ++xp) xsum = fma(sx + xp >= 0 && sx + xp < src_cols ? - CONVERT_TO_WT(loadpix(srcptr + mad24(xp, TSIZE, src_index))) : scalar, coeffs_x[xp], xsum); + CONVERT_TO_WT(loadpix(srcptr + mad24(xp, TSIZE, src_index))) : scalar, (WT)(coeffs_x[xp]), xsum); } - sum = fma(xsum, coeffs_y[yp], sum); + sum = fma(xsum, (WT)(coeffs_y[yp]), sum); } else - sum = fma(scalar, coeffs_y[yp], sum); + sum = fma(scalar, (WT)(coeffs_y[yp]), sum); } storepix(CONVERT_TO_T(sum), dst); diff --git a/modules/imgproc/test/ocl/test_warp.cpp b/modules/imgproc/test/ocl/test_warp.cpp index d44ffd102e..eac9548ea0 100644 --- a/modules/imgproc/test/ocl/test_warp.cpp +++ b/modules/imgproc/test/ocl/test_warp.cpp @@ -52,8 +52,12 @@ //M*/ #include "../test_precomp.hpp" +#include "opencv2/core/types.hpp" #include "opencv2/ts/ocl_test.hpp" +#include +#include + #ifdef HAVE_OPENCL namespace opencv_test { @@ -458,6 +462,279 @@ OCL_TEST_P(Remap_INTER_LINEAR, Mat) } } +template +class RemapBorderConstant : public ::testing::Test +{ +public: + RemapBorderConstant() + { + constexpr float k1 = -0.08f; + constexpr float coordinate_offset = -0.03f; + + undistortedCanvas = calculateUndistortedCanvas(source_size, k1); + mapx.create(undistortedCanvas.size()); + mapy.create(undistortedCanvas.size()); + + for (int y = 0; y < undistortedCanvas.height; ++y) + { + for (int x = 0; x < undistortedCanvas.width; ++x) + { + const Point2f undistorted(undistortedCanvas.tl() + cv::Point2i(x, y)); + const Point2f distorted = barrelDistort(undistorted, source_size, + k1); + mapx.at(y, x) = x == 0 ? coordinate_offset : distorted.x; + mapy.at(y, x) = y == 0 ? coordinate_offset : distorted.y; + } + } + } + +protected: + void runTensorProductTest() + { + const Mat_ source = createTensorProductImage(source_size); + runTest(createBarrelDistortedImage(source, barrel_coefficient)); + } + + void runCheckerboardTest() + { + const Mat_ source = createCheckerboardImage(source_size); + runTest(createBarrelDistortedImage(source, barrel_coefficient)); + } + +private: + static const Size source_size; + static constexpr float barrel_coefficient = -0.08f; + + void runTest(const Mat_& src) + { + constexpr float border_value = 2.0f; + constexpr float epsilon = 0.02f; + const Mat_ expected = remapReference(src, mapx, mapy, + border_value); + + Mat_ cpuResult; + OCL_OFF(cv::remap(src, cpuResult, mapx, mapy, INTER_LINEAR, + BORDER_CONSTANT, Scalar::all(border_value))); + + UMat usrc, umapx, umapy, uResult; + OCL_ON(src.copyTo(usrc)); + OCL_ON(mapx.copyTo(umapx)); + OCL_ON(mapy.copyTo(umapy)); + OCL_ON(cv::remap(usrc, uResult, umapx, umapy, INTER_LINEAR, + BORDER_CONSTANT, Scalar::all(border_value))); + Mat_ oclResult; + uResult.copyTo(oclResult); + + EXPECT_MAT_NEAR(expected, cpuResult, epsilon); + EXPECT_MAT_NEAR(expected, oclResult, epsilon); + EXPECT_MAT_NEAR(cpuResult, oclResult, epsilon); + } + + static Mat_ createTensorProductImage(const Size& image_size) + { + constexpr float first_coordinate = 0.0f; + constexpr float last_coordinate = 1.0f; + Mat_ image(image_size); + + for (int y = 0; y < image.rows; ++y) + { + const T yValue = first_coordinate + + static_cast(y) * (last_coordinate - first_coordinate) / + static_cast(image.rows - 1); + + for (int x = 0; x < image.cols; ++x) + { + const T xValue = first_coordinate + + static_cast(x) * (last_coordinate - first_coordinate) / + static_cast(image.cols - 1); + image.template at(y, x) = xValue * yValue; + } + } + + return image; + } + + static Mat_ createCheckerboardImage(const Size& image_size) + { + constexpr int checkerboard_period = 2; + constexpr float first_coordinate = 0.0f; + constexpr float last_coordinate = 1.0f; + + Mat_ image(image_size); + + for (int y = 0; y < image.rows; ++y) + { + for (int x = 0; x < image.cols; ++x) + { + image.template at(y, x) = (x + y) % checkerboard_period == 0 ? + first_coordinate : last_coordinate; + } + } + + return image; + } + + static Mat_ remapReference(const Mat_& src, const Mat1f& mapx, + const Mat1f& mapy, T borderValue) + { + CV_Assert(src.type() == DataType::type); + CV_Assert(mapx.size() == mapy.size()); + + Mat_ dst(mapx.size()); + + for (int y = 0; y < dst.rows; ++y) + { + for (int x = 0; x < dst.cols; ++x) + { + const float fx = mapx.at(y, x); + const float fy = mapy.at(y, x); + const int scaledX = static_cast(std::floor(fx * INTER_TAB_SIZE + 0.5f)); + const int scaledY = static_cast(std::floor(fy * INTER_TAB_SIZE + 0.5f)); + const int sx = scaledX >> INTER_BITS; + const int sy = scaledY >> INTER_BITS; + const T wx = static_cast(scaledX & (INTER_TAB_SIZE - 1)) / + static_cast(INTER_TAB_SIZE); + const T wy = static_cast(scaledY & (INTER_TAB_SIZE - 1)) / + static_cast(INTER_TAB_SIZE); + + T value = 0; + + for (int yp = 0; yp < 2; ++yp) + { + for (int xp = 0; xp < 2; ++xp) + { + const int sourceX = sx + xp; + const int sourceY = sy + yp; + const bool inBounds = sourceX >= 0 && sourceX < src.cols && + sourceY >= 0 && sourceY < src.rows; + const T sourceValue = inBounds ? src.template at(sourceY, sourceX) : borderValue; + const T xWeight = xp == 0 ? 1 - wx : wx; + const T yWeight = yp == 0 ? 1 - wy : wy; + value += sourceValue * xWeight * yWeight; + } + } + + dst.template at(y, x) = value; + } + } + + return dst; + } + + static Point2f barrelDistort(const Point2f& point, const Size& imageSize, float k1) + { + const Point2f center = cv::Point2i(imageSize - cv::Size(1, 1)) / 2.0f; + const float scale = std::max(center.x, center.y); + const Point2f normalized = (point - center) / scale; + const float radiusSquared = normalized.dot(normalized); + const float factor = 1.0f + k1 * radiusSquared; + + return center + normalized * scale * factor; + } + + static Point2f barrelUndistort(const Point2f& point, const Size& imageSize, float k1) + { + const Point2f center = cv::Point2i(imageSize - cv::Size(1, 1)) / 2.0f; + const float scale = std::max(center.x, center.y); + const Point2f normalized = (point - center) / scale; + const float distortedRadius = std::hypot(normalized.x, normalized.y); + + if (std::fpclassify(distortedRadius) == FP_ZERO) + return point; + + float undistortedRadius = distortedRadius; + constexpr int inverse_iterations = 8; + + for (int i = 0; i < inverse_iterations; ++i) + { + const float radiusSquared = undistortedRadius * undistortedRadius; + const float function = undistortedRadius * (1.0f + k1 * radiusSquared) - + distortedRadius; + const float derivative = 1.0f + 3.0f * k1 * radiusSquared; + undistortedRadius -= function / derivative; + } + + const float radialScale = undistortedRadius / distortedRadius; + return center + normalized * radialScale * scale; + } + + static Rect calculateUndistortedCanvas(const Size& sourceSize, float k1) + { + const Point2f corners[] = { + Point2f(0.0f, 0.0f), + Point2f(static_cast(sourceSize.width - 1), 0.0f), + Point2f(0.0f, static_cast(sourceSize.height - 1)), + Point2f(static_cast(sourceSize.width - 1), + static_cast(sourceSize.height - 1)) + }; + + float minX = std::numeric_limits::max(); + float minY = std::numeric_limits::max(); + float maxX = std::numeric_limits::lowest(); + float maxY = std::numeric_limits::lowest(); + + for (const Point2f& corner : corners) + { + const Point2f undistortedCorner = barrelUndistort(corner, sourceSize, k1); + minX = std::min(minX, undistortedCorner.x); + minY = std::min(minY, undistortedCorner.y); + maxX = std::max(maxX, undistortedCorner.x); + maxY = std::max(maxY, undistortedCorner.y); + } + + const int left = static_cast(std::floor(minX)); + const int top = static_cast(std::floor(minY)); + const int right = static_cast(std::ceil(maxX)); + const int bottom = static_cast(std::ceil(maxY)); + + return Rect(left, top, right - left + 1, bottom - top + 1); + } + + static Mat_ createBarrelDistortedImage(const Mat_& image, float k1) + { + Mat1f mapx(image.size()); + Mat1f mapy(image.size()); + + for (int y = 0; y < image.rows; ++y) + { + for (int x = 0; x < image.cols; ++x) + { + const Point2f undistorted = barrelUndistort(Point2i(x, y), + image.size(), k1); + mapx.at(y, x) = undistorted.x; + mapy.at(y, x) = undistorted.y; + } + } + + constexpr float first_coordinate = 0.0f; + return remapReference(image, mapx, mapy, first_coordinate); + } + + Mat1f mapx; + Mat1f mapy; + Rect undistortedCanvas; +}; + +template +const Size RemapBorderConstant::source_size{30, 20}; +template +constexpr float RemapBorderConstant::barrel_coefficient; + +using RemapBorderConstantTypes = ::testing::Types; +TYPED_TEST_CASE(RemapBorderConstant, RemapBorderConstantTypes); + +// The smooth tensor-product image exposes artifacts at the image edges. +TYPED_TEST(RemapBorderConstant, tensor_product) +{ + this->runTensorProductTest(); +} + +// The checkerboard image exposes interpolation artifacts within the image. +TYPED_TEST(RemapBorderConstant, checkerboard) +{ + this->runCheckerboardTest(); +} + ///////////////////////////////////////////////////////////////////////////////////////////////// // remap relative