From 4fe51e51e075f1c0e1efbc1b7b2ff92a4d222013 Mon Sep 17 00:00:00 2001 From: Madan mohan Manokar Date: Tue, 14 Jul 2026 17:47:22 +0530 Subject: [PATCH] Merge pull request #29507 from amd:fast_remap_ext imgproc: Optimized remap interpolation #29507 - Add a SIMD dispatch file for remap and vectorize the single-channel (C1) in-bounds paths of bilinear, bicubic and lanczos4 interpolation (32F / 16U / 16S) using width-agnostic gather - Dispatch the bilinear C1 path and drop the per-pixel weight-table gathers - Widen the fixed-point coordinate map conversion - 32F lanczos4 is kept on the scalar path: its vectorized 64-tap accumulation deviates beyond the set float accuracy tolerance - Vectorize the bilinear inlier/outlier run detection so any-channel linear remap skips constant-status runs with SIMD instead of a per-pixel bounds test ### 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 --- modules/imgproc/CMakeLists.txt | 1 + modules/imgproc/src/imgwarp.avx2.cpp | 4 + modules/imgproc/src/imgwarp.cpp | 259 +++++++++++--- modules/imgproc/src/imgwarp.simd.hpp | 478 +++++++++++++++++++++++++ modules/imgproc/src/imgwarp.sse4_1.cpp | 4 + 5 files changed, 706 insertions(+), 40 deletions(-) create mode 100644 modules/imgproc/src/imgwarp.simd.hpp diff --git a/modules/imgproc/CMakeLists.txt b/modules/imgproc/CMakeLists.txt index 6ff2621c80..8fac3e83bc 100644 --- a/modules/imgproc/CMakeLists.txt +++ b/modules/imgproc/CMakeLists.txt @@ -11,6 +11,7 @@ ocv_add_dispatched_file(morph SSE2 SSE4_1 AVX2) ocv_add_dispatched_file(smooth SSE2 SSE4_1 AVX2 AVX512_ICL) ocv_add_dispatched_file(sumpixels SSE2 AVX2 AVX512_SKX) ocv_add_dispatched_file(equalize_hist AVX512_ICL) +ocv_add_dispatched_file(imgwarp SSE4_1 AVX2 AVX512_SKX AVX512_ICL) ocv_define_module(imgproc opencv_core WRAP java objc python js) if(OPENCV_CORE_EXCLUDE_C_API) diff --git a/modules/imgproc/src/imgwarp.avx2.cpp b/modules/imgproc/src/imgwarp.avx2.cpp index 6a795a7c92..393ac707a2 100644 --- a/modules/imgproc/src/imgwarp.avx2.cpp +++ b/modules/imgproc/src/imgwarp.avx2.cpp @@ -13,6 +13,7 @@ // Copyright (C) 2000-2008, Intel Corporation, all rights reserved. // Copyright (C) 2009, Willow Garage Inc., all rights reserved. // Copyright (C) 2014-2015, Itseez Inc., all rights reserved. +// Copyright (C) 2026, Advanced Micro Devices, all rights reserved. // Third party copyrights are property of their respective owners. // // Redistribution and use in source and binary forms, with or without modification, @@ -95,4 +96,7 @@ int warpAffineBlockline(int *adelta, int *bdelta, short* xy, short* alpha, int X } } + +#include "imgwarp.simd.hpp" + /* End of file. */ diff --git a/modules/imgproc/src/imgwarp.cpp b/modules/imgproc/src/imgwarp.cpp index ee5b62fda9..17f77c1599 100644 --- a/modules/imgproc/src/imgwarp.cpp +++ b/modules/imgproc/src/imgwarp.cpp @@ -13,6 +13,7 @@ // Copyright (C) 2000-2008, Intel Corporation, all rights reserved. // Copyright (C) 2009, Willow Garage Inc., all rights reserved. // Copyright (C) 2014-2015, Itseez Inc., all rights reserved. +// Copyright (C) 2026, Advanced Micro Devices, all rights reserved. // Third party copyrights are property of their respective owners. // // Redistribution and use in source and binary forms, with or without modification, @@ -55,6 +56,9 @@ #include "opencv2/core/softfloat.hpp" #include "imgwarp.hpp" +#include "imgwarp.simd.hpp" +#include "imgwarp.simd_declarations.hpp" // defines CV_CPU_DISPATCH_MODES_ALL=AVX512_ICL,...,BASELINE based on CMakeLists.txt content + using namespace cv; namespace cv @@ -611,6 +615,81 @@ template using RemapVec_8u = RemapNoVec; #endif +template +struct RemapBilinearVecC1 +{ + int operator()(const T*, size_t, T*, const short*, const ushort*, + const AT*, int, int, int) const { return 0; } +}; + +template<> +struct RemapBilinearVecC1 +{ + int operator()(const float* S0, size_t sstep, float* D, const short* XY, + const ushort* FXY, const float* wtab, int dx, int X1, int off_y) const + { + CV_CPU_DISPATCH(remapBilinearC1_simd, + (CV_32F, (const uchar*)S0, sstep, (uchar*)D, XY, FXY, wtab, dx, X1, off_y), + CV_CPU_DISPATCH_MODES_ALL); + } +}; + +template<> +struct RemapBilinearVecC1 +{ + int operator()(const ushort* S0, size_t sstep, ushort* D, const short* XY, + const ushort* FXY, const float* wtab, int dx, int X1, int off_y) const + { + CV_CPU_DISPATCH(remapBilinearC1_simd, + (CV_16U, (const uchar*)S0, sstep, (uchar*)D, XY, FXY, wtab, dx, X1, off_y), + CV_CPU_DISPATCH_MODES_ALL); + } +}; + +template<> +struct RemapBilinearVecC1 +{ + int operator()(const short* S0, size_t sstep, short* D, const short* XY, + const ushort* FXY, const float* wtab, int dx, int X1, int off_y) const + { + CV_CPU_DISPATCH(remapBilinearC1_simd, + (CV_16S, (const uchar*)S0, sstep, (uchar*)D, XY, FXY, wtab, dx, X1, off_y), + CV_CPU_DISPATCH_MODES_ALL); + } +}; + +static inline int remapBilinearSameRun( const short* XY, int dx, int end, + unsigned width1, unsigned height1, bool inl ) +{ + int n = 0; +#if (CV_SIMD || CV_SIMD_SCALABLE) + const int span = VTraits::vlanes(); + const v_int16 vw = vx_setall_s16((short)std::min(width1, 0x7fff)); + const v_int16 vh = vx_setall_s16((short)std::min(height1, 0x7fff)); + const v_int16 vm1 = vx_setall_s16(-1); + for( ; dx + n + span <= end; n += span ) + { + v_int16 sx, sy; + v_load_deinterleave(XY + (dx + n) * 2, sx, sy); + // in-bounds: 0 <= sx < width1 && 0 <= sy < height1 + v_int16 inb = v_and(v_and(v_gt(sx, vm1), v_lt(sx, vw)), + v_and(v_gt(sy, vm1), v_lt(sy, vh))); + const bool allSame = inl ? v_check_all(inb) : !v_check_any(inb); + if( !allSame ) + break; + } + vx_cleanup(); +#endif + for( ; dx + n < end; n++ ) + { + const int sx = XY[(dx + n) * 2], sy = XY[(dx + n) * 2 + 1]; + const bool ib = (unsigned)sx < width1 && (unsigned)sy < height1; + if( ib != inl ) + break; + } + return n; +} + template static void remapBilinear( const Mat& _src, Mat& _dst, const Mat& _xy, const Mat& _fxy, const void* _wtab, @@ -647,6 +726,12 @@ static void remapBilinear( const Mat& _src, Mat& _dst, const Mat& _xy, const int off_y = (isRelative ? (_offset.y+dy) : 0); for(int dx = 0; dx <= dsize.width; dx++ ) { + if( !isRelative && dx < dsize.width ) + { + int n = remapBilinearSameRun(XY, dx, dsize.width, width1, height1, prevInlier); + if( n > 0 ) + dx += n - 1; + } bool curInlier = dx < dsize.width ? (unsigned)XY[dx*2]+(isRelative ? (_offset.x+dx) : 0) < width1 && (unsigned)XY[dx*2+1]+off_y < height1 : !prevInlier; @@ -667,6 +752,11 @@ static void remapBilinear( const Mat& _src, Mat& _dst, const Mat& _xy, if( cn == 1 ) { + if( !isRelative ) + { + int n = RemapBilinearVecC1()(S0, sstep, D, XY, FXY, wtab, dx, X1, off_y); + D += n; dx += n; + } for( ; dx < X1; dx++, D++ ) { int sx = XY[dx*2]+(isRelative ? (_offset.x+dx) : 0), sy = XY[dx*2+1]+off_y; @@ -843,6 +933,53 @@ static void remapBilinear( const Mat& _src, Mat& _dst, const Mat& _xy, } +// Dispatch shim for the single-channel non-relative bicubic in-bounds fast path (32F only). +template +struct RemapBicubicVecC1 +{ + int operator()(const T*, size_t, T*, const short*, const ushort*, const AT*, + int, int, unsigned, unsigned, int) const { return 0; } +}; + +template<> +struct RemapBicubicVecC1 +{ + int operator()(const float* S0, size_t sstep, float* D, const short* XY, + const ushort* FXY, const float* wtab, int dx, int dwidth, unsigned width1, + unsigned height1, int off_y) const + { + CV_CPU_DISPATCH(remapBicubicC1wp_simd, + (CV_32F, (const uchar*)S0, sstep, (uchar*)D, XY, FXY, wtab, dx, dwidth, width1, height1, off_y), + CV_CPU_DISPATCH_MODES_ALL); + } +}; + +template<> +struct RemapBicubicVecC1 +{ + int operator()(const ushort* S0, size_t sstep, ushort* D, const short* XY, + const ushort* FXY, const float* wtab, int dx, int dwidth, unsigned width1, + unsigned height1, int off_y) const + { + CV_CPU_DISPATCH(remapBicubicC1wp_simd, + (CV_16U, (const uchar*)S0, sstep, (uchar*)D, XY, FXY, wtab, dx, dwidth, width1, height1, off_y), + CV_CPU_DISPATCH_MODES_ALL); + } +}; + +template<> +struct RemapBicubicVecC1 +{ + int operator()(const short* S0, size_t sstep, short* D, const short* XY, + const ushort* FXY, const float* wtab, int dx, int dwidth, unsigned width1, + unsigned height1, int off_y) const + { + CV_CPU_DISPATCH(remapBicubicC1wp_simd, + (CV_16S, (const uchar*)S0, sstep, (uchar*)D, XY, FXY, wtab, dx, dwidth, width1, height1, off_y), + CV_CPU_DISPATCH_MODES_ALL); + } +}; + template static void remapBicubic( const Mat& _src, Mat& _dst, const Mat& _xy, const Mat& _fxy, const void* _wtab, @@ -879,6 +1016,12 @@ static void remapBicubic( const Mat& _src, Mat& _dst, const Mat& _xy, const int off_y = isRelative ? (_offset.y+dy) : 0; for(int dx = 0; dx < dsize.width; dx++, D += cn ) { + if( cn == 1 && !isRelative ) + { + int n = RemapBicubicVecC1()(S0, sstep, D, XY, FXY, wtab, dx, + dsize.width, width1, height1, off_y); + if( n > 0 ) { D += (n - 1)*cn; dx += n - 1; continue; } + } 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; @@ -948,6 +1091,33 @@ static void remapBicubic( const Mat& _src, Mat& _dst, const Mat& _xy, } +template +struct RemapLanczos4VecC1 +{ + int operator()(const T*, size_t, T*, const short*, const ushort*, const AT*, + int, int, unsigned, unsigned, int) const { return 0; } +}; + +#define CV_REMAP_LANCZOS4_SHIM(T, DEPTH) \ +template<> struct RemapLanczos4VecC1 \ +{ \ + int operator()(const T* S0, size_t sstep, T* D, const short* XY, \ + const ushort* FXY, const float* wtab, int dx, int dwidth, \ + unsigned width1, unsigned height1, int off_y) const \ + { \ + CV_CPU_DISPATCH(remapLanczos4C1_simd, \ + (DEPTH, (const uchar*)S0, sstep, (uchar*)D, XY, FXY, wtab, \ + dx, dwidth, width1, height1, off_y), \ + CV_CPU_DISPATCH_MODES_ALL); \ + } \ +}; +// 32F is intentionally not shimmed: its vectorized accumulation deviates beyond +// the float accuracy tolerance, so it stays on the scalar loop. Emitting a shim +// would add a per-pixel dispatch call that returns 0 and slows the scalar path. +CV_REMAP_LANCZOS4_SHIM(ushort, CV_16U) +CV_REMAP_LANCZOS4_SHIM(short, CV_16S) +#undef CV_REMAP_LANCZOS4_SHIM + template static void remapLanczos4( const Mat& _src, Mat& _dst, const Mat& _xy, const Mat& _fxy, const void* _wtab, @@ -984,6 +1154,12 @@ static void remapLanczos4( const Mat& _src, Mat& _dst, const Mat& _xy, const int off_y = isRelative ? (_offset.y+dy) : 0; for(int dx = 0; dx < dsize.width; dx++, D += cn ) { + if( cn == 1 && !isRelative ) + { + int n = RemapLanczos4VecC1()(S0, sstep, D, XY, FXY, wtab, dx, + dsize.width, width1, height1, off_y); + if( n > 0 ) { D += (n - 1)*cn; dx += n - 1; continue; } + } const int off_x = isRelative ? (_offset.x+dx) : 0; int sx = XY[dx*2]-3+off_x, sy = XY[dx*2+1]-3+off_y; const AT* w = wtab + FXY[dx]*64; @@ -1131,21 +1307,21 @@ public: const float* sY = m2->ptr(y+y1) + x; x1 = 0; - #if CV_SIMD128 + #if (CV_SIMD || CV_SIMD_SCALABLE) { - int span = VTraits::vlanes(); + int span = VTraits::vlanes(); for( ; x1 <= bcols - span * 2; x1 += span * 2 ) { - v_int32x4 ix0 = v_round(v_load(sX + x1)); - v_int32x4 iy0 = v_round(v_load(sY + x1)); - v_int32x4 ix1 = v_round(v_load(sX + x1 + span)); - v_int32x4 iy1 = v_round(v_load(sY + x1 + span)); + v_int32 ix0 = v_round(vx_load(sX + x1)); + v_int32 iy0 = v_round(vx_load(sY + x1)); + v_int32 ix1 = v_round(vx_load(sX + x1 + span)); + v_int32 iy1 = v_round(vx_load(sY + x1 + span)); - v_int16x8 dx, dy; - dx = v_pack(ix0, ix1); - dy = v_pack(iy0, iy1); + v_int16 dx = v_pack(ix0, ix1); + v_int16 dy = v_pack(iy0, iy1); v_store_interleave(XY + x1 * 2, dx, dy); } + vx_cleanup(); } #endif for( ; x1 < bcols; x1++ ) @@ -1172,12 +1348,13 @@ public: const ushort* sA = m2->ptr(y+y1) + x; x1 = 0; - #if CV_SIMD128 + #if (CV_SIMD || CV_SIMD_SCALABLE) { - v_uint16x8 v_scale = v_setall_u16(INTER_TAB_SIZE2 - 1); - int span = VTraits::vlanes(); + v_uint16 v_scale = vx_setall_u16(INTER_TAB_SIZE2 - 1); + int span = VTraits::vlanes(); for( ; x1 <= bcols - span; x1 += span ) - v_store((unsigned short*)(A + x1), v_and(v_load(sA + x1), v_scale)); + v_store((unsigned short*)(A + x1), v_and(vx_load(sA + x1), v_scale)); + vx_cleanup(); } #endif for( ; x1 < bcols; x1++ ) @@ -1189,26 +1366,27 @@ public: const float* sY = m2->ptr(y+y1) + x; x1 = 0; - #if CV_SIMD128 + #if (CV_SIMD || CV_SIMD_SCALABLE) { - v_float32x4 v_scale = v_setall_f32((float)INTER_TAB_SIZE); - v_int32x4 v_scale2 = v_setall_s32(INTER_TAB_SIZE - 1); - int span = VTraits::vlanes(); + v_float32 v_scale = vx_setall_f32((float)INTER_TAB_SIZE); + v_int32 v_scale2 = vx_setall_s32(INTER_TAB_SIZE - 1); + int span = VTraits::vlanes(); for( ; x1 <= bcols - span * 2; x1 += span * 2 ) { - v_int32x4 v_sx0 = v_round(v_mul(v_scale, v_load(sX + x1))); - v_int32x4 v_sy0 = v_round(v_mul(v_scale, v_load(sY + x1))); - v_int32x4 v_sx1 = v_round(v_mul(v_scale, v_load(sX + x1 + span))); - v_int32x4 v_sy1 = v_round(v_mul(v_scale, v_load(sY + x1 + span))); - v_uint16x8 v_sx8 = v_reinterpret_as_u16(v_pack(v_and(v_sx0, v_scale2), v_and(v_sx1, v_scale2))); - v_uint16x8 v_sy8 = v_reinterpret_as_u16(v_pack(v_and(v_sy0, v_scale2), v_and(v_sy1, v_scale2))); - v_uint16x8 v_v = v_or(v_shl(v_sy8), v_sx8); + v_int32 v_sx0 = v_round(v_mul(v_scale, vx_load(sX + x1))); + v_int32 v_sy0 = v_round(v_mul(v_scale, vx_load(sY + x1))); + v_int32 v_sx1 = v_round(v_mul(v_scale, vx_load(sX + x1 + span))); + v_int32 v_sy1 = v_round(v_mul(v_scale, vx_load(sY + x1 + span))); + v_uint16 v_sx8 = v_reinterpret_as_u16(v_pack(v_and(v_sx0, v_scale2), v_and(v_sx1, v_scale2))); + v_uint16 v_sy8 = v_reinterpret_as_u16(v_pack(v_and(v_sy0, v_scale2), v_and(v_sy1, v_scale2))); + v_uint16 v_v = v_or(v_shl(v_sy8), v_sx8); v_store(A + x1, v_v); - v_int16x8 v_d0 = v_pack(v_shr(v_sx0), v_shr(v_sx1)); - v_int16x8 v_d1 = v_pack(v_shr(v_sy0), v_shr(v_sy1)); + v_int16 v_d0 = v_pack(v_shr(v_sx0), v_shr(v_sx1)); + v_int16 v_d1 = v_pack(v_shr(v_sy0), v_shr(v_sy1)); v_store_interleave(XY + (x1 << 1), v_d0, v_d1); } + vx_cleanup(); } #endif for( ; x1 < bcols; x1++ ) @@ -1226,28 +1404,29 @@ public: const float* sXY = m1->ptr(y+y1) + x*2; x1 = 0; - #if CV_SIMD128 + #if (CV_SIMD || CV_SIMD_SCALABLE) { - v_float32x4 v_scale = v_setall_f32((float)INTER_TAB_SIZE); - v_int32x4 v_scale2 = v_setall_s32(INTER_TAB_SIZE - 1), v_scale3 = v_setall_s32(INTER_TAB_SIZE); - int span = VTraits::vlanes(); + v_float32 v_scale = vx_setall_f32((float)INTER_TAB_SIZE); + v_int32 v_scale2 = vx_setall_s32(INTER_TAB_SIZE - 1), v_scale3 = vx_setall_s32(INTER_TAB_SIZE); + int span = VTraits::vlanes(); for( ; x1 <= bcols - span * 2; x1 += span * 2 ) { - v_float32x4 v_fx, v_fy; + v_float32 v_fx, v_fy; v_load_deinterleave(sXY + (x1 << 1), v_fx, v_fy); - v_int32x4 v_sx0 = v_round(v_mul(v_fx, v_scale)); - v_int32x4 v_sy0 = v_round(v_mul(v_fy, v_scale)); + v_int32 v_sx0 = v_round(v_mul(v_fx, v_scale)); + v_int32 v_sy0 = v_round(v_mul(v_fy, v_scale)); v_load_deinterleave(sXY + ((x1 + span) << 1), v_fx, v_fy); - v_int32x4 v_sx1 = v_round(v_mul(v_fx, v_scale)); - v_int32x4 v_sy1 = v_round(v_mul(v_fy, v_scale)); - v_int32x4 v_v0 = v_muladd(v_scale3, (v_and(v_sy0, v_scale2)), (v_and(v_sx0, v_scale2))); - v_int32x4 v_v1 = v_muladd(v_scale3, (v_and(v_sy1, v_scale2)), (v_and(v_sx1, v_scale2))); - v_uint16x8 v_v8 = v_reinterpret_as_u16(v_pack(v_v0, v_v1)); + v_int32 v_sx1 = v_round(v_mul(v_fx, v_scale)); + v_int32 v_sy1 = v_round(v_mul(v_fy, v_scale)); + v_int32 v_v0 = v_muladd(v_scale3, (v_and(v_sy0, v_scale2)), (v_and(v_sx0, v_scale2))); + v_int32 v_v1 = v_muladd(v_scale3, (v_and(v_sy1, v_scale2)), (v_and(v_sx1, v_scale2))); + v_uint16 v_v8 = v_reinterpret_as_u16(v_pack(v_v0, v_v1)); v_store(A + x1, v_v8); - v_int16x8 v_dx = v_pack(v_shr(v_sx0), v_shr(v_sx1)); - v_int16x8 v_dy = v_pack(v_shr(v_sy0), v_shr(v_sy1)); + v_int16 v_dx = v_pack(v_shr(v_sx0), v_shr(v_sx1)); + v_int16 v_dy = v_pack(v_shr(v_sy0), v_shr(v_sy1)); v_store_interleave(XY + (x1 << 1), v_dx, v_dy); } + vx_cleanup(); } #endif diff --git a/modules/imgproc/src/imgwarp.simd.hpp b/modules/imgproc/src/imgwarp.simd.hpp new file mode 100644 index 0000000000..f7eb4397b7 --- /dev/null +++ b/modules/imgproc/src/imgwarp.simd.hpp @@ -0,0 +1,478 @@ +// 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. +// +// Copyright (C) 2026, Advanced Micro Devices, all rights reserved. + +#include "opencv2/core/hal/intrin.hpp" + +namespace cv { + +CV_CPU_OPTIMIZATION_NAMESPACE_BEGIN + +int remapBilinearC1_simd(int depth, const uchar* S0, size_t sstep, uchar* D, + const short* XY, const ushort* FXY, const float* wtab, + int dx, int X1, int off_y); + +int remapBicubicC1_simd(int depth, const uchar* S0, size_t sstep, uchar* D, + const short* XY, const ushort* FXY, + int dx, int dwidth, unsigned width1, unsigned height1, int off_y); + +int remapLanczos4C1_simd(int depth, const uchar* S0, size_t sstep, uchar* D, + const short* XY, const ushort* FXY, const float* wtab, + int dx, int dwidth, unsigned width1, unsigned height1, int off_y); + + +int remapBicubicC1wp_simd(int depth, const uchar* S0, size_t sstep, uchar* D, + const short* XY, const ushort* FXY, const float* wtab, + int dx, int dwidth, unsigned width1, unsigned height1, int off_y); + +#ifndef CV_CPU_OPTIMIZATION_DECLARATIONS_ONLY + +#if (CV_SIMD || CV_SIMD_SCALABLE) + +static inline v_float32 remapGatherF32(const float* base, const int* ofs) +{ + float CV_DECL_ALIGNED(CV_SIMD_WIDTH) buf[VTraits::max_nlanes]; + const int n = VTraits::vlanes(); + for (int k = 0; k < n; k++) + buf[k] = base[ofs[k]]; + return vx_load(buf); +} + +static inline v_float32 remapGatherF32(const ushort* base, const int* ofs) +{ + float CV_DECL_ALIGNED(CV_SIMD_WIDTH) buf[VTraits::max_nlanes]; + const int n = VTraits::vlanes(); + for (int k = 0; k < n; k++) + buf[k] = (float)base[ofs[k]]; + return vx_load(buf); +} + +static inline v_float32 remapGatherF32(const short* base, const int* ofs) +{ + float CV_DECL_ALIGNED(CV_SIMD_WIDTH) buf[VTraits::max_nlanes]; + const int n = VTraits::vlanes(); + for (int k = 0; k < n; k++) + buf[k] = (float)base[ofs[k]]; + return vx_load(buf); +} + +static CV_ALWAYS_INLINE void remapCorners(const short* S0, size_t sstep, const int* ofs, + v_float32& s0, v_float32& s1, + v_float32& s2, v_float32& s3) +{ + int CV_DECL_ALIGNED(CV_SIMD_WIDTH) topbuf[VTraits::max_nlanes]; + int CV_DECL_ALIGNED(CV_SIMD_WIDTH) botbuf[VTraits::max_nlanes]; + const int n = VTraits::vlanes(); + for (int k = 0; k < n; k++) + { + const short* p = S0 + ofs[k]; + int t, b; + memcpy(&t, p, sizeof(t)); + memcpy(&b, p + sstep, sizeof(b)); + topbuf[k] = t; botbuf[k] = b; + } + v_int32 top = vx_load(topbuf), bot = vx_load(botbuf); + s0 = v_cvt_f32(v_shr<16>(v_shl<16>(top))); // low 16 bits (sign-extended) + s1 = v_cvt_f32(v_shr<16>(top)); // high 16 bits (sign-extended) + s2 = v_cvt_f32(v_shr<16>(v_shl<16>(bot))); + s3 = v_cvt_f32(v_shr<16>(bot)); +} + +static CV_ALWAYS_INLINE void remapCorners(const ushort* S0, size_t sstep, const int* ofs, + v_float32& s0, v_float32& s1, + v_float32& s2, v_float32& s3) +{ + int CV_DECL_ALIGNED(CV_SIMD_WIDTH) topbuf[VTraits::max_nlanes]; + int CV_DECL_ALIGNED(CV_SIMD_WIDTH) botbuf[VTraits::max_nlanes]; + const int n = VTraits::vlanes(); + for (int k = 0; k < n; k++) + { + const ushort* p = S0 + ofs[k]; + int t, b; + memcpy(&t, p, sizeof(t)); + memcpy(&b, p + sstep, sizeof(b)); + topbuf[k] = t; botbuf[k] = b; + } + const v_uint32 lo16 = vx_setall_u32(0xffff); + v_uint32 top = v_reinterpret_as_u32(vx_load(topbuf)); + v_uint32 bot = v_reinterpret_as_u32(vx_load(botbuf)); + s0 = v_cvt_f32(v_reinterpret_as_s32(v_and(top, lo16))); // low 16 (zero-ext) + s1 = v_cvt_f32(v_reinterpret_as_s32(v_shr<16>(top))); // high 16 (zero-ext) + s2 = v_cvt_f32(v_reinterpret_as_s32(v_and(bot, lo16))); + s3 = v_cvt_f32(v_reinterpret_as_s32(v_shr<16>(bot))); +} + +static inline void remapStoreC1(float* D, const v_float32& res) +{ + v_store(D, res); +} + +static inline void remapStoreC1(ushort* D, const v_float32& res) +{ + v_pack_u_store(D, v_round(res)); +} + +static inline void remapStoreC1(short* D, const v_float32& res) +{ + v_pack_store(D, v_round(res)); +} + +template +static int remapBilinearC1_run(const T* S0, size_t sstep, T* D, + const short* XY, const ushort* FXY, + const float* wtab, int dx, int X1, int off_y) +{ + CV_UNUSED(wtab); + const int vlanes = VTraits::vlanes(); + const int dx0 = dx; + const v_float32 vone = vx_setall_f32(1.f); + const v_float32 vscale = vx_setall_f32(1.f / INTER_TAB_SIZE); + const v_int32 vmask = vx_setall_s32(INTER_TAB_SIZE - 1); + int CV_DECL_ALIGNED(CV_SIMD_WIDTH) ofs[VTraits::max_nlanes]; + for( ; dx <= X1 - vlanes; dx += vlanes ) + { + for( int k = 0; k < vlanes; k++ ) + { + const int sx = XY[(dx + k) * 2]; + const int sy = XY[(dx + k) * 2 + 1] + off_y; + ofs[k] = sy * (int)sstep + sx; + } + v_float32 s0, s1, s2, s3; + remapCorners(S0, sstep, ofs, s0, s1, s2, s3); + + v_int32 fxy = v_reinterpret_as_s32(vx_load_expand(FXY + dx)); + v_float32 fx = v_mul(v_cvt_f32(v_and(fxy, vmask)), vscale); + v_float32 fy = v_mul(v_cvt_f32(v_shr(fxy)), vscale); + v_float32 cx0 = v_sub(vone, fx); + v_float32 cy0 = v_sub(vone, fy); + v_float32 w0 = v_mul(cx0, cy0); + v_float32 w1 = v_mul(fx, cy0); + v_float32 w2 = v_mul(cx0, fy); + v_float32 w3 = v_mul(fx, fy); + + v_float32 res = v_fma(s0, w0, v_fma(s1, w1, v_fma(s2, w2, v_mul(s3, w3)))); + remapStoreC1(D + (dx - dx0), res); + } + vx_cleanup(); + return dx - dx0; +} + +static int remapBilinearF32_run(const float* S0, size_t sstep, float* D, + const short* XY, const ushort* FXY, + const float* wtab, int dx, int X1, int off_y) +{ + CV_UNUSED(wtab); + const int vlanes = VTraits::vlanes(); + const int dx0 = dx; + const float* S1 = S0 + sstep; + const v_float32 vone = vx_setall_f32(1.f); + const v_float32 vscale = vx_setall_f32(1.f / INTER_TAB_SIZE); + const v_int32 vmask = vx_setall_s32(INTER_TAB_SIZE - 1); + int CV_DECL_ALIGNED(CV_SIMD_WIDTH) ofs[VTraits::max_nlanes]; + for( ; dx <= X1 - vlanes; dx += vlanes ) + { + for( int k = 0; k < vlanes; k++ ) + { + const int sx = XY[(dx + k) * 2]; + const int sy = XY[(dx + k) * 2 + 1] + off_y; + ofs[k] = sy * (int)sstep + sx; + } + v_float32 s0 = remapGatherF32(S0, ofs); + v_float32 s1 = remapGatherF32(S0 + 1, ofs); + v_float32 s2 = remapGatherF32(S1, ofs); + v_float32 s3 = remapGatherF32(S1 + 1, ofs); + + v_int32 fxy = v_reinterpret_as_s32(vx_load_expand(FXY + dx)); + v_float32 fx = v_mul(v_cvt_f32(v_and(fxy, vmask)), vscale); + v_float32 fy = v_mul(v_cvt_f32(v_shr(fxy)), vscale); + v_float32 cx0 = v_sub(vone, fx); + v_float32 cy0 = v_sub(vone, fy); + v_float32 w0 = v_mul(cx0, cy0); + v_float32 w1 = v_mul(fx, cy0); + v_float32 w2 = v_mul(cx0, fy); + v_float32 w3 = v_mul(fx, fy); + + v_float32 res = v_fma(s0, w0, v_fma(s1, w1, v_fma(s2, w2, v_mul(s3, w3)))); + v_store(D + (dx - dx0), res); + } + vx_cleanup(); + return dx - dx0; +} + +// Evaluate the four cubic interpolation coefficients for a vector of fractional +// positions, matching interpolateCubic() (A = -0.75). The strict remap test +// tolerates an absolute error of 1.0 for bicubic, so FMA contraction is fine. +static inline void interpolateCubicV(const v_float32& x, + v_float32& c0, v_float32& c1, + v_float32& c2, v_float32& c3) +{ + const v_float32 A = vx_setall_f32(-0.75f); + const v_float32 A5 = vx_setall_f32(-3.75f); // 5*A + const v_float32 A8 = vx_setall_f32(-6.0f); // 8*A + const v_float32 A4 = vx_setall_f32(-3.0f); // 4*A + const v_float32 Ap2 = vx_setall_f32(1.25f); // A+2 + const v_float32 Ap3 = vx_setall_f32(2.25f); // A+3 + const v_float32 one = vx_setall_f32(1.f); + + v_float32 xp1 = v_add(x, one); + c0 = v_sub(v_mul(v_add(v_mul(v_sub(v_mul(A, xp1), A5), xp1), A8), xp1), A4); + c1 = v_add(v_mul(v_mul(v_sub(v_mul(Ap2, x), Ap3), x), x), one); + v_float32 u = v_sub(one, x); + c2 = v_add(v_mul(v_mul(v_sub(v_mul(Ap2, u), Ap3), u), u), one); + c3 = v_sub(v_sub(v_sub(one, c0), c1), c2); +} + +// Select one of four vectors by runtime index. Used instead of an array of +// vector types, which is not valid for sizeless RVV vector types. +static inline v_float32 selectV4(int i, const v_float32& a, const v_float32& b, + const v_float32& c, const v_float32& d) +{ + return i == 0 ? a : (i == 1 ? b : (i == 2 ? c : d)); +} + +template +static int remapBicubicC1_run(const T* S0, size_t sstep, T* D, const short* XY, + const ushort* FXY, int dx, int dwidth, + unsigned width1, unsigned height1, int off_y) +{ + const int vlanes = VTraits::vlanes(); + const int dx0 = dx; + const v_float32 vscale = vx_setall_f32(1.f / INTER_TAB_SIZE); + const v_int32 vmask = vx_setall_s32(INTER_TAB_SIZE - 1); + int CV_DECL_ALIGNED(CV_SIMD_WIDTH) ofs[VTraits::max_nlanes]; + float CV_DECL_ALIGNED(CV_SIMD_WIDTH) buf[VTraits::max_nlanes]; + for( ; dx <= dwidth - vlanes; dx += vlanes ) + { + bool allIn = true; + for( int k = 0; k < vlanes; k++ ) + { + const unsigned sx = (unsigned)(XY[(dx + k) * 2] - 1); + const unsigned sy = (unsigned)(XY[(dx + k) * 2 + 1] - 1 + off_y); + if( sx >= width1 || sy >= height1 ) { allIn = false; break; } + ofs[k] = (int)((XY[(dx + k) * 2 + 1] - 1 + off_y) * (int)sstep + + (XY[(dx + k) * 2] - 1)); + } + if( !allIn ) + break; + + v_int32 fxy = v_reinterpret_as_s32(vx_load_expand(FXY + dx)); + v_float32 fx = v_mul(v_cvt_f32(v_and(fxy, vmask)), vscale); + v_float32 fy = v_mul(v_cvt_f32(v_shr(fxy)), vscale); + v_float32 vx0, vx1, vx2, vx3, vy0, vy1, vy2, vy3; + interpolateCubicV(fx, vx0, vx1, vx2, vx3); + interpolateCubicV(fy, vy0, vy1, vy2, vy3); + + v_float32 acc = vx_setzero_f32(); + for( int r = 0; r < 4; r++ ) + { + const int roff = r * (int)sstep; + const v_float32 vyr = selectV4(r, vy0, vy1, vy2, vy3); + for( int c = 0; c < 4; c++ ) + { + for( int k = 0; k < vlanes; k++ ) + buf[k] = (float)S0[ofs[k] + roff + c]; + acc = v_fma(vx_load(buf), v_mul(vyr, selectV4(c, vx0, vx1, vx2, vx3)), acc); + } + } + remapStoreC1(D + (dx - dx0), acc); + } + vx_cleanup(); + return dx - dx0; +} + +#if CV_SIMD128 +static inline void remapLoad8(const float* S, v_float32x4& lo, v_float32x4& hi) +{ + lo = v_load(S); hi = v_load(S + 4); +} +static inline void remapLoad8(const ushort* S, v_float32x4& lo, v_float32x4& hi) +{ + v_uint16x8 v = v_load(S); + v_uint32x4 a, b; v_expand(v, a, b); + lo = v_cvt_f32(v_reinterpret_as_s32(a)); + hi = v_cvt_f32(v_reinterpret_as_s32(b)); +} +static inline void remapLoad8(const short* S, v_float32x4& lo, v_float32x4& hi) +{ + v_int16x8 v = v_load(S); + v_int32x4 a, b; v_expand(v, a, b); + lo = v_cvt_f32(a); hi = v_cvt_f32(b); +} + +static inline void remapStoreScalar(float* D, float v) { *D = v; } +static inline void remapStoreScalar(ushort* D, float v) { *D = saturate_cast(v); } +static inline void remapStoreScalar(short* D, float v) { *D = saturate_cast(v); } + +static inline v_float32x4 remapLoad4(const float* S) { return v_load(S); } +static inline v_float32x4 remapLoad4(const ushort* S) { return v_cvt_f32(v_reinterpret_as_s32(v_load_expand(S))); } +static inline v_float32x4 remapLoad4(const short* S) { return v_cvt_f32(v_load_expand(S)); } + +template +static int remapBicubicC1wp_run(const T* S0, size_t sstep, T* D, const short* XY, + const ushort* FXY, const float* wtab, int dx, + int dwidth, unsigned width1, unsigned height1, int off_y) +{ + const int dx0 = dx; + for( ; dx < dwidth; dx++ ) + { + const unsigned sx = (unsigned)(XY[dx * 2] - 1); + const unsigned sy = (unsigned)(XY[dx * 2 + 1] - 1 + off_y); + if( sx >= width1 || sy >= height1 ) + break; + const float* w = wtab + FXY[dx] * 16; + const T* S = S0 + (size_t)sy * sstep + sx; + v_float32x4 acc = v_setzero_f32(); + for( int r = 0; r < 4; r++, S += sstep, w += 4 ) + acc = v_fma(remapLoad4(S), v_load(w), acc); + remapStoreScalar(D + (dx - dx0), v_reduce_sum(acc)); + } + vx_cleanup(); + return dx - dx0; +} + +template +static int remapLanczos4C1_run(const T* S0, size_t sstep, T* D, const short* XY, + const ushort* FXY, const float* wtab, int dx, + int dwidth, unsigned width1, unsigned height1, int off_y) +{ + const int dx0 = dx; + for( ; dx < dwidth; dx++ ) + { + const unsigned sx = (unsigned)(XY[dx * 2] - 3); + const unsigned sy = (unsigned)(XY[dx * 2 + 1] - 3 + off_y); + if( sx >= width1 || sy >= height1 ) + break; + const float* w = wtab + FXY[dx] * 64; + const T* S = S0 + (size_t)sy * sstep + sx; + v_float32x4 acc = v_setzero_f32(); + for( int r = 0; r < 8; r++, S += sstep, w += 8 ) + { + v_float32x4 s_lo, s_hi; + remapLoad8(S, s_lo, s_hi); + acc = v_fma(s_lo, v_load(w), acc); + acc = v_fma(s_hi, v_load(w + 4), acc); + } + remapStoreScalar(D + (dx - dx0), v_reduce_sum(acc)); + } + vx_cleanup(); + return dx - dx0; +} +#endif // CV_SIMD128 + +#endif // CV_SIMD + +int remapBilinearC1_simd(int depth, const uchar* S0, size_t sstep, uchar* D, + const short* XY, const ushort* FXY, const float* wtab, + int dx, int X1, int off_y) +{ +#if (CV_SIMD || CV_SIMD_SCALABLE) + switch (depth) + { + case CV_32F: + return remapBilinearF32_run((const float*)S0, sstep, (float*)D, + XY, FXY, wtab, dx, X1, off_y); + case CV_16U: + return remapBilinearC1_run((const ushort*)S0, sstep, (ushort*)D, + XY, FXY, wtab, dx, X1, off_y); + case CV_16S: + return remapBilinearC1_run((const short*)S0, sstep, (short*)D, + XY, FXY, wtab, dx, X1, off_y); + default: + return 0; + } +#else + CV_UNUSED(depth); CV_UNUSED(S0); CV_UNUSED(sstep); CV_UNUSED(D); + CV_UNUSED(XY); CV_UNUSED(FXY); CV_UNUSED(wtab); + CV_UNUSED(dx); CV_UNUSED(X1); CV_UNUSED(off_y); + return 0; +#endif +} + +int remapBicubicC1_simd(int depth, const uchar* S0, size_t sstep, uchar* D, + const short* XY, const ushort* FXY, + int dx, int dwidth, unsigned width1, unsigned height1, int off_y) +{ +#if (CV_SIMD || CV_SIMD_SCALABLE) + switch (depth) + { + case CV_32F: + return remapBicubicC1_run((const float*)S0, sstep, (float*)D, + XY, FXY, dx, dwidth, width1, height1, off_y); + case CV_16U: + return remapBicubicC1_run((const ushort*)S0, sstep, (ushort*)D, + XY, FXY, dx, dwidth, width1, height1, off_y); + case CV_16S: + return remapBicubicC1_run((const short*)S0, sstep, (short*)D, + XY, FXY, dx, dwidth, width1, height1, off_y); + default: + return 0; + } +#else + CV_UNUSED(depth); CV_UNUSED(S0); CV_UNUSED(sstep); CV_UNUSED(D); + CV_UNUSED(XY); CV_UNUSED(FXY); CV_UNUSED(dx); CV_UNUSED(dwidth); + CV_UNUSED(width1); CV_UNUSED(height1); CV_UNUSED(off_y); + return 0; +#endif +} + +int remapLanczos4C1_simd(int depth, const uchar* S0, size_t sstep, uchar* D, + const short* XY, const ushort* FXY, const float* wtab, + int dx, int dwidth, unsigned width1, unsigned height1, int off_y) +{ +#if CV_SIMD128 + switch (depth) + { + // CV_32F is intentionally omitted: the vectorized 8x8 accumulation reorders + // the 64-tap sum, so on the tight 1e-3 float tolerance it diverges from the + // scalar path (used for relative maps). 32F lanczos4 stays on the scalar loop. + case CV_16U: + return remapLanczos4C1_run((const ushort*)S0, sstep, (ushort*)D, + XY, FXY, wtab, dx, dwidth, width1, height1, off_y); + case CV_16S: + return remapLanczos4C1_run((const short*)S0, sstep, (short*)D, + XY, FXY, wtab, dx, dwidth, width1, height1, off_y); + default: + return 0; + } +#else + CV_UNUSED(depth); CV_UNUSED(S0); CV_UNUSED(sstep); CV_UNUSED(D); + CV_UNUSED(XY); CV_UNUSED(FXY); CV_UNUSED(wtab); CV_UNUSED(dx); CV_UNUSED(dwidth); + CV_UNUSED(width1); CV_UNUSED(height1); CV_UNUSED(off_y); + return 0; +#endif +} + +int remapBicubicC1wp_simd(int depth, const uchar* S0, size_t sstep, uchar* D, + const short* XY, const ushort* FXY, const float* wtab, + int dx, int dwidth, unsigned width1, unsigned height1, int off_y) +{ +#if CV_SIMD128 + switch (depth) + { + case CV_32F: + return remapBicubicC1wp_run((const float*)S0, sstep, (float*)D, + XY, FXY, wtab, dx, dwidth, width1, height1, off_y); + case CV_16U: + return remapBicubicC1wp_run((const ushort*)S0, sstep, (ushort*)D, + XY, FXY, wtab, dx, dwidth, width1, height1, off_y); + case CV_16S: + return remapBicubicC1wp_run((const short*)S0, sstep, (short*)D, + XY, FXY, wtab, dx, dwidth, width1, height1, off_y); + default: + return 0; + } +#else + CV_UNUSED(depth); CV_UNUSED(S0); CV_UNUSED(sstep); CV_UNUSED(D); + CV_UNUSED(XY); CV_UNUSED(FXY); CV_UNUSED(wtab); CV_UNUSED(dx); CV_UNUSED(dwidth); + CV_UNUSED(width1); CV_UNUSED(height1); CV_UNUSED(off_y); + return 0; +#endif +} + +#endif // CV_CPU_OPTIMIZATION_DECLARATIONS_ONLY + +CV_CPU_OPTIMIZATION_NAMESPACE_END + +} // namespace cv diff --git a/modules/imgproc/src/imgwarp.sse4_1.cpp b/modules/imgproc/src/imgwarp.sse4_1.cpp index 9f7ac82fca..1038d17ab5 100644 --- a/modules/imgproc/src/imgwarp.sse4_1.cpp +++ b/modules/imgproc/src/imgwarp.sse4_1.cpp @@ -13,6 +13,7 @@ // Copyright (C) 2000-2008, Intel Corporation, all rights reserved. // Copyright (C) 2009, Willow Garage Inc., all rights reserved. // Copyright (C) 2014-2015, Itseez Inc., all rights reserved. +// Copyright (C) 2026, Advanced Micro Devices, all rights reserved. // Third party copyrights are property of their respective owners. // // Redistribution and use in source and binary forms, with or without modification, @@ -466,4 +467,7 @@ Ptr WarpPerspectiveLine_SSE4::getImpl(const double *M) } } + +#include "imgwarp.simd.hpp" + /* End of file. */