diff --git a/modules/core/CMakeLists.txt b/modules/core/CMakeLists.txt index 414dc9907b..f5ebf3e720 100644 --- a/modules/core/CMakeLists.txt +++ b/modules/core/CMakeLists.txt @@ -8,7 +8,7 @@ ocv_add_dispatched_file(convert_scale SSE2 AVX2 LASX) ocv_add_dispatched_file(count_non_zero SSE2 AVX2 LASX) ocv_add_dispatched_file(has_non_zero SSE2 AVX2 LASX ) ocv_add_dispatched_file(matmul SSE2 SSE4_1 AVX2 AVX512_SKX NEON_DOTPROD LASX) -ocv_add_dispatched_file(mean SSE2 AVX2 LASX) +ocv_add_dispatched_file(mean SSE2 AVX2 AVX512_SKX AVX512_ICL LASX) ocv_add_dispatched_file(merge SSE2 AVX2 LASX) ocv_add_dispatched_file(split SSE2 AVX2 LASX) ocv_add_dispatched_file(sum SSE2 AVX2 LASX) diff --git a/modules/core/src/mean.simd.hpp b/modules/core/src/mean.simd.hpp index fa0ac7d3cd..972fe9a24d 100644 --- a/modules/core/src/mean.simd.hpp +++ b/modules/core/src/mean.simd.hpp @@ -2,6 +2,7 @@ // 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, Inc., all rights reserved. #include "precomp.hpp" #include "stat.hpp" @@ -31,7 +32,35 @@ struct SumSqr_SIMD { int operator () (const uchar * src0, const uchar * mask, int * sum, int * sqsum, int len, int cn) const { - if (mask || (cn != 1 && cn != 2 && cn != 4)) + if (mask) + return 0; + // cn==3: deinterleave so each lane belongs to one channel, then reduce per + // channel. sum/sqsum stay exact in s32 within the dispatcher's 1<<15 block. + if (cn == 3) + { + const int vl8 = VTraits::vlanes(); + v_int32 vs0 = vx_setzero_s32(), vs1 = vx_setzero_s32(), vs2 = vx_setzero_s32(); + v_int32 vq0 = vx_setzero_s32(), vq1 = vx_setzero_s32(), vq2 = vx_setzero_s32(); + auto acc = [](const v_uint8& ch, v_int32& vs, v_int32& vq) + { + v_uint16 lo, hi; v_expand(ch, lo, hi); + v_uint32 s0, s1, s2, s3; v_expand(lo, s0, s1); v_expand(hi, s2, s3); + vs = v_add(vs, v_reinterpret_as_s32(v_add(v_add(s0, s1), v_add(s2, s3)))); + v_int16 d0 = v_reinterpret_as_s16(lo), d1 = v_reinterpret_as_s16(hi); + vq = v_add(vq, v_add(v_dotprod(d0, d0), v_dotprod(d1, d1))); + }; + int e = 0; + for (; e <= len - vl8; e += vl8) + { + v_uint8 a, b, c; v_load_deinterleave(src0 + e * 3, a, b, c); + acc(a, vs0, vq0); acc(b, vs1, vq1); acc(c, vs2, vq2); + } + sum[0] += v_reduce_sum(vs0); sum[1] += v_reduce_sum(vs1); sum[2] += v_reduce_sum(vs2); + sqsum[0] += v_reduce_sum(vq0); sqsum[1] += v_reduce_sum(vq1); sqsum[2] += v_reduce_sum(vq2); + v_cleanup(); + return e; + } + if (cn != 1 && cn != 2 && cn != 4) return 0; len *= cn; @@ -98,7 +127,32 @@ struct SumSqr_SIMD { int operator () (const schar * src0, const uchar * mask, int * sum, int * sqsum, int len, int cn) const { - if (mask || (cn != 1 && cn != 2 && cn != 4)) + if (mask) + return 0; + if (cn == 3) + { + const int vl8 = VTraits::vlanes(); + v_int32 vs0 = vx_setzero_s32(), vs1 = vx_setzero_s32(), vs2 = vx_setzero_s32(); + v_int32 vq0 = vx_setzero_s32(), vq1 = vx_setzero_s32(), vq2 = vx_setzero_s32(); + auto acc = [](const v_int8& ch, v_int32& vs, v_int32& vq) + { + v_int16 lo, hi; v_expand(ch, lo, hi); + v_int32 s0, s1, s2, s3; v_expand(lo, s0, s1); v_expand(hi, s2, s3); + vs = v_add(vs, v_add(v_add(s0, s1), v_add(s2, s3))); + vq = v_add(vq, v_add(v_dotprod(lo, lo), v_dotprod(hi, hi))); + }; + int e = 0; + for (; e <= len - vl8; e += vl8) + { + v_int8 a, b, c; v_load_deinterleave(src0 + e * 3, a, b, c); + acc(a, vs0, vq0); acc(b, vs1, vq1); acc(c, vs2, vq2); + } + sum[0] += v_reduce_sum(vs0); sum[1] += v_reduce_sum(vs1); sum[2] += v_reduce_sum(vs2); + sqsum[0] += v_reduce_sum(vq0); sqsum[1] += v_reduce_sum(vq1); sqsum[2] += v_reduce_sum(vq2); + v_cleanup(); + return e; + } + if (cn != 1 && cn != 2 && cn != 4) return 0; len *= cn; @@ -160,6 +214,212 @@ struct SumSqr_SIMD } }; +#if (CV_SIMD_64F || CV_SIMD_SCALABLE_64F) + +// Expand a 16-bit vector into two s32 halves. The unsigned overload reinterprets +// after a widening expand (values 0..65535 stay non-negative in s32, so the +// subsequent v_cvt_f64 is exact); the signed overload expands directly. This is +// the only real difference between the ushort and short reductions below. +static inline void sumSqrExpandS32(const v_uint16& v, v_int32& lo, v_int32& hi) +{ + v_uint32 a, b; v_expand(v, a, b); + lo = v_reinterpret_as_s32(a); hi = v_reinterpret_as_s32(b); +} +static inline void sumSqrExpandS32(const v_int16& v, v_int32& lo, v_int32& hi) +{ + v_expand(v, lo, hi); +} + +// Shared 16-bit -> {s32 sum, f64 sqsum} reduction for ushort/short. The register +// type is deduced from the pointer, so the body is identical for both depths. +template +static inline int sumSqr16ToF64(const T* src0, int* sum, double* sqsum, int len, int cn) +{ + typedef decltype(vx_load(src0)) VT; // v_uint16 or v_int16 + + if (cn == 3) + { + const int vl = VTraits::vlanes(); + v_int32 vs0 = vx_setzero_s32(), vs1 = vx_setzero_s32(), vs2 = vx_setzero_s32(); + v_float64 vq0 = vx_setzero_f64(), vq1 = vx_setzero_f64(), vq2 = vx_setzero_f64(); + auto acc = [](const VT& ch, v_int32& vs, v_float64& vq) + { + v_int32 lo, hi; sumSqrExpandS32(ch, lo, hi); + vs = v_add(vs, v_add(lo, hi)); + v_float64 f0 = v_cvt_f64(lo), f1 = v_cvt_f64_high(lo); + v_float64 f2 = v_cvt_f64(hi), f3 = v_cvt_f64_high(hi); + vq = v_fma(f0, f0, v_fma(f1, f1, v_fma(f2, f2, v_fma(f3, f3, vq)))); + }; + int e = 0; + for (; e <= len - vl; e += vl) + { + VT a, b, c; v_load_deinterleave(src0 + e * 3, a, b, c); + acc(a, vs0, vq0); acc(b, vs1, vq1); acc(c, vs2, vq2); + } + sum[0] += v_reduce_sum(vs0); sum[1] += v_reduce_sum(vs1); sum[2] += v_reduce_sum(vs2); + sqsum[0] += v_reduce_sum(vq0); sqsum[1] += v_reduce_sum(vq1); sqsum[2] += v_reduce_sum(vq2); + v_cleanup(); + return e; + } + if (cn != 1 && cn != 2 && cn != 4) + return 0; + len *= cn; + + const int vl16 = VTraits::vlanes(); + const int vl32 = VTraits::vlanes(); + const int vl64 = VTraits::vlanes(); + + // The lane->channel scatter below (i % cn) and the returned pixel count + // (x / cn) require the lane counts to be a multiple of cn. This always holds + // for fixed-width SIMD (vlanes is a power of two >= cn for cn in {1,2,4}), + // but scalable backends (e.g. RVV) may report a non-multiple; fall back to + // scalar there. vl16 == 2*vl32, so checking vl32 covers both scatter loops. + if (vl32 % cn != 0) + return 0; + + v_int32 v_sum = vx_setzero_s32(); + v_float64 v_sq0 = vx_setzero_f64(), v_sq1 = vx_setzero_f64(); + v_float64 v_sq2 = vx_setzero_f64(), v_sq3 = vx_setzero_f64(); + + int x = 0; + for (; x <= len - vl16; x += vl16) + { + VT v_src = vx_load(src0 + x); + v_int32 lo, hi; sumSqrExpandS32(v_src, lo, hi); + v_sum = v_add(v_sum, v_add(lo, hi)); + + v_float64 f0 = v_cvt_f64(lo); + v_float64 f1 = v_cvt_f64_high(lo); + v_float64 f2 = v_cvt_f64(hi); + v_float64 f3 = v_cvt_f64_high(hi); + v_sq0 = v_fma(f0, f0, v_sq0); + v_sq1 = v_fma(f1, f1, v_sq1); + v_sq2 = v_fma(f2, f2, v_sq2); + v_sq3 = v_fma(f3, f3, v_sq3); + } + + int CV_DECL_ALIGNED(CV_SIMD_WIDTH) sbuf[VTraits::max_nlanes]; + double CV_DECL_ALIGNED(CV_SIMD_WIDTH) qbuf[VTraits::max_nlanes * 4]; + v_store_aligned(sbuf, v_sum); + v_store_aligned(qbuf + vl64 * 0, v_sq0); + v_store_aligned(qbuf + vl64 * 1, v_sq1); + v_store_aligned(qbuf + vl64 * 2, v_sq2); + v_store_aligned(qbuf + vl64 * 3, v_sq3); + + // sbuf lane j (j channel j%cn (vl32%cn==0). + // qbuf index i (i channel i%cn (vl16%cn==0). + for (int i = 0; i < vl32; ++i) + sum[i % cn] += sbuf[i]; + for (int i = 0; i < vl16; ++i) + sqsum[i % cn] += qbuf[i]; + + v_cleanup(); + return x / cn; +} + +// Shared 32-bit -> f64 reduction for int/float. Both sum and sqsum accumulate in +// f64 to match the scalar reference's double accumulation; values are widened to +// f64 before squaring so the square is computed in double (same as (double)v*v). +// The register type is deduced from the pointer, so int and float share one body. +template +static inline int sumSqr32ToF64(const T* src0, double* sum, double* sqsum, int len, int cn) +{ + typedef decltype(vx_load(src0)) VT; // v_int32 or v_float32 + + if (cn == 3) + { + const int vl = VTraits::vlanes(); + v_float64 vs0 = vx_setzero_f64(), vs1 = vx_setzero_f64(), vs2 = vx_setzero_f64(); + v_float64 vq0 = vx_setzero_f64(), vq1 = vx_setzero_f64(), vq2 = vx_setzero_f64(); + auto acc = [](const VT& ch, v_float64& vs, v_float64& vq) + { + v_float64 f0 = v_cvt_f64(ch), f1 = v_cvt_f64_high(ch); + vs = v_add(vs, v_add(f0, f1)); + vq = v_fma(f0, f0, v_fma(f1, f1, vq)); + }; + int e = 0; + for (; e <= len - vl; e += vl) + { + VT a, b, c; v_load_deinterleave(src0 + e * 3, a, b, c); + acc(a, vs0, vq0); acc(b, vs1, vq1); acc(c, vs2, vq2); + } + sum[0] += v_reduce_sum(vs0); sum[1] += v_reduce_sum(vs1); sum[2] += v_reduce_sum(vs2); + sqsum[0] += v_reduce_sum(vq0); sqsum[1] += v_reduce_sum(vq1); sqsum[2] += v_reduce_sum(vq2); + v_cleanup(); + return e; + } + if (cn != 1 && cn != 2 && cn != 4) + return 0; + len *= cn; + + const int vl32 = VTraits::vlanes(); + const int vl64 = VTraits::vlanes(); + + // See note in sumSqr16ToF64: guard scalable backends whose lane count may + // not divide cn (the i % cn scatter / x / cn rely on it). + if (vl32 % cn != 0) + return 0; + + v_float64 vs0 = vx_setzero_f64(), vs1 = vx_setzero_f64(); + v_float64 vq0 = vx_setzero_f64(), vq1 = vx_setzero_f64(); + + int x = 0; + for (; x <= len - vl32; x += vl32) + { + VT v_src = vx_load(src0 + x); + v_float64 f0 = v_cvt_f64(v_src); + v_float64 f1 = v_cvt_f64_high(v_src); + vs0 = v_add(vs0, f0); + vs1 = v_add(vs1, f1); + vq0 = v_fma(f0, f0, vq0); + vq1 = v_fma(f1, f1, vq1); + } + + double CV_DECL_ALIGNED(CV_SIMD_WIDTH) sbuf[VTraits::max_nlanes * 2]; + double CV_DECL_ALIGNED(CV_SIMD_WIDTH) qbuf[VTraits::max_nlanes * 2]; + v_store_aligned(sbuf, vs0); v_store_aligned(sbuf + vl64, vs1); + v_store_aligned(qbuf, vq0); v_store_aligned(qbuf + vl64, vq1); + + for (int i = 0; i < vl32; ++i) + { + sum[i % cn] += sbuf[i]; + sqsum[i % cn] += qbuf[i]; + } + + v_cleanup(); + return x / cn; +} + +template <> +struct SumSqr_SIMD +{ + int operator () (const ushort* src0, const uchar* mask, int* sum, double* sqsum, int len, int cn) const + { return mask ? 0 : sumSqr16ToF64(src0, sum, sqsum, len, cn); } +}; + +template <> +struct SumSqr_SIMD +{ + int operator () (const short* src0, const uchar* mask, int* sum, double* sqsum, int len, int cn) const + { return mask ? 0 : sumSqr16ToF64(src0, sum, sqsum, len, cn); } +}; + +template <> +struct SumSqr_SIMD +{ + int operator () (const int* src0, const uchar* mask, double* sum, double* sqsum, int len, int cn) const + { return mask ? 0 : sumSqr32ToF64(src0, sum, sqsum, len, cn); } +}; + +template <> +struct SumSqr_SIMD +{ + int operator () (const float* src0, const uchar* mask, double* sum, double* sqsum, int len, int cn) const + { return mask ? 0 : sumSqr32ToF64(src0, sum, sqsum, len, cn); } +}; + +#endif // CV_SIMD_64F + #endif template @@ -253,6 +513,21 @@ static int sumsqr_(const T* src0, const uchar* mask, ST* sum, SQT* sqsum, int le sum[0] = s0; sqsum[0] = sq0; } + else if( cn == 2 ) + { + ST s0 = sum[0], s1 = sum[1]; + SQT sq0 = sqsum[0], sq1 = sqsum[1]; + for( i = 0; i < len; i++, src += 2 ) + if( mask[i] ) + { + T v0 = src[0], v1 = src[1]; + s0 += v0; sq0 += (SQT)v0*v0; + s1 += v1; sq1 += (SQT)v1*v1; + nzm++; + } + sum[0] = s0; sum[1] = s1; + sqsum[0] = sq0; sqsum[1] = sq1; + } else if( cn == 3 ) { ST s0 = sum[0], s1 = sum[1], s2 = sum[2]; @@ -269,6 +544,23 @@ static int sumsqr_(const T* src0, const uchar* mask, ST* sum, SQT* sqsum, int le sum[0] = s0; sum[1] = s1; sum[2] = s2; sqsum[0] = sq0; sqsum[1] = sq1; sqsum[2] = sq2; } + else if( cn == 4 ) + { + ST s0 = sum[0], s1 = sum[1], s2 = sum[2], s3 = sum[3]; + SQT sq0 = sqsum[0], sq1 = sqsum[1], sq2 = sqsum[2], sq3 = sqsum[3]; + for( i = 0; i < len; i++, src += 4 ) + if( mask[i] ) + { + T v0 = src[0], v1 = src[1], v2 = src[2], v3 = src[3]; + s0 += v0; sq0 += (SQT)v0*v0; + s1 += v1; sq1 += (SQT)v1*v1; + s2 += v2; sq2 += (SQT)v2*v2; + s3 += v3; sq3 += (SQT)v3*v3; + nzm++; + } + sum[0] = s0; sum[1] = s1; sum[2] = s2; sum[3] = s3; + sqsum[0] = sq0; sqsum[1] = sq1; sqsum[2] = sq2; sqsum[3] = sq3; + } else { for( i = 0; i < len; i++, src += cn )