diff --git a/modules/core/src/mathfuncs.cpp b/modules/core/src/mathfuncs.cpp index c3bd50e0d4..8060f0bd33 100644 --- a/modules/core/src/mathfuncs.cpp +++ b/modules/core/src/mathfuncs.cpp @@ -1360,20 +1360,32 @@ bool checkRange(InputArray _src, bool quiet, Point* pt, double minVal, double ma #ifdef HAVE_OPENCL -static bool ocl_patchNaNs( InputOutputArray _a, float value ) +static bool ocl_patchNaNs( InputOutputArray _a, double value ) { + int ftype = _a.depth(); + + bool doubleSupport = ocl::Device::getDefault().doubleFPConfig() > 0; + if (!doubleSupport && ftype == CV_64F) + return false; + int rowsPerWI = ocl::Device::getDefault().isIntel() ? 4 : 1; ocl::Kernel k("KF", ocl::core::arithm_oclsrc, - format("-D UNARY_OP -D OP_PATCH_NANS -D dstT=float -D DEPTH_dst=%d -D rowsPerWI=%d", - CV_32F, rowsPerWI)); + format("-D UNARY_OP -D OP_PATCH_NANS -D dstT=%s -D DEPTH_dst=%d -D rowsPerWI=%d %s", + ftype == CV_64F ? "double" : "float", ftype, rowsPerWI, + doubleSupport ? "-D DOUBLE_SUPPORT" : "")); if (k.empty()) return false; UMat a = _a.getUMat(); int cn = a.channels(); - k.args(ocl::KernelArg::ReadOnlyNoSize(a), - ocl::KernelArg::WriteOnly(a, cn), (float)value); + // the value is passed with the same type as the kernel operates on + if (ftype == CV_64F) + k.args(ocl::KernelArg::ReadOnlyNoSize(a), + ocl::KernelArg::WriteOnly(a, cn), value); + else + k.args(ocl::KernelArg::ReadOnlyNoSize(a), + ocl::KernelArg::WriteOnly(a, cn), (float)value); size_t globalsize[2] = { (size_t)a.cols * cn, ((size_t)a.rows + rowsPerWI - 1) / rowsPerWI }; return k.run(2, globalsize, NULL, false); @@ -1381,55 +1393,99 @@ static bool ocl_patchNaNs( InputOutputArray _a, float value ) #endif +static void patchNaNs_32f( uchar* ptr, size_t ulen, double newVal ) +{ + int* tptr = (int*)ptr; + int len = (int)ulen; + int j = 0; + Cv32suf val; + val.f = (float)newVal; + +#if (CV_SIMD || CV_SIMD_SCALABLE) + v_int32 v_pos_mask = vx_setall_s32(0x7fffffff), v_exp_mask = vx_setall_s32(0x7f800000); + v_int32 v_val = vx_setall_s32(val.i); + + int cWidth = VTraits::vlanes(); + for (; j < len - cWidth * 2 + 1; j += cWidth * 2) + { + v_int32 v_src0 = vx_load(tptr + j); + v_int32 v_src1 = vx_load(tptr + j + cWidth); + + v_int32 v_cmp_mask0 = v_lt(v_exp_mask, v_and(v_src0, v_pos_mask)); + v_int32 v_cmp_mask1 = v_lt(v_exp_mask, v_and(v_src1, v_pos_mask)); + + if (v_check_any(v_or(v_cmp_mask0, v_cmp_mask1))) + { + v_int32 v_dst0 = v_select(v_cmp_mask0, v_val, v_src0); + v_int32 v_dst1 = v_select(v_cmp_mask1, v_val, v_src1); + + v_store(tptr + j, v_dst0); + v_store(tptr + j + cWidth, v_dst1); + } + } +#endif + + for( ; j < len; j++ ) + if( (tptr[j] & 0x7fffffff) > 0x7f800000 ) + tptr[j] = val.i; +} + +static void patchNaNs_64f( uchar* ptr, size_t ulen, double newVal ) +{ + double* tptr = (double*)ptr; + int len = (int)ulen; + int j = 0; + + // universal intrinsics have no ordering comparison for 64-bit integers, + // so NaNs are detected with v_not_nan() rather than with the + // exponent/mantissa trick used for CV_32F above +#if (CV_SIMD_64F || CV_SIMD_SCALABLE_64F) + v_float64 v_val = vx_setall_f64(newVal); + + int cWidth = VTraits::vlanes(); + for (; j < len - cWidth * 2 + 1; j += cWidth * 2) + { + v_float64 v_src0 = vx_load(tptr + j); + v_float64 v_src1 = vx_load(tptr + j + cWidth); + + v_float64 v_ok_mask0 = v_not_nan(v_src0); + v_float64 v_ok_mask1 = v_not_nan(v_src1); + + if (!v_check_all(v_ok_mask0) || !v_check_all(v_ok_mask1)) + { + v_store(tptr + j, v_select(v_ok_mask0, v_src0, v_val)); + v_store(tptr + j + cWidth, v_select(v_ok_mask1, v_src1, v_val)); + } + } +#endif + + for( ; j < len; j++ ) + if( cvIsNaN(tptr[j]) ) + tptr[j] = newVal; +} + void patchNaNs( InputOutputArray _a, double _val ) { CV_INSTRUMENT_REGION(); - CV_Assert( _a.depth() == CV_32F ); + int depth = _a.depth(); + CV_Assert( depth == CV_32F || depth == CV_64F ); CV_OCL_RUN(_a.isUMat() && _a.dims() <= 2, - ocl_patchNaNs(_a, (float)_val)) + ocl_patchNaNs(_a, _val)) Mat a = _a.getMat(); const Mat* arrays[] = {&a, 0}; - int* ptrs[1] = {}; - NAryMatIterator it(arrays, (uchar**)ptrs); - int len = (int)(it.size*a.channels()); - Cv32suf val; - val.f = (float)_val; + uchar* ptrs[1] = {}; + NAryMatIterator it(arrays, ptrs); + size_t len = it.size*a.channels(); for( size_t i = 0; i < it.nplanes; i++, ++it ) { - int* tptr = ptrs[0]; - int j = 0; - -#if (CV_SIMD || CV_SIMD_SCALABLE) - v_int32 v_pos_mask = vx_setall_s32(0x7fffffff), v_exp_mask = vx_setall_s32(0x7f800000); - v_int32 v_val = vx_setall_s32(val.i); - - int cWidth = VTraits::vlanes(); - for (; j < len - cWidth * 2 + 1; j += cWidth * 2) - { - v_int32 v_src0 = vx_load(tptr + j); - v_int32 v_src1 = vx_load(tptr + j + cWidth); - - v_int32 v_cmp_mask0 = v_lt(v_exp_mask, v_and(v_src0, v_pos_mask)); - v_int32 v_cmp_mask1 = v_lt(v_exp_mask, v_and(v_src1, v_pos_mask)); - - if (v_check_any(v_or(v_cmp_mask0, v_cmp_mask1))) - { - v_int32 v_dst0 = v_select(v_cmp_mask0, v_val, v_src0); - v_int32 v_dst1 = v_select(v_cmp_mask1, v_val, v_src1); - - v_store(tptr + j, v_dst0); - v_store(tptr + j + cWidth, v_dst1); - } - } -#endif - - for( ; j < len; j++ ) - if( (tptr[j] & 0x7fffffff) > 0x7f800000 ) - tptr[j] = val.i; + if( depth == CV_64F ) + patchNaNs_64f(ptrs[0], len, _val); + else + patchNaNs_32f(ptrs[0], len, _val); } } diff --git a/modules/core/test/ocl/test_arithm.cpp b/modules/core/test/ocl/test_arithm.cpp index 2f20091c32..3067b3c444 100644 --- a/modules/core/test/ocl/test_arithm.cpp +++ b/modules/core/test/ocl/test_arithm.cpp @@ -1721,8 +1721,35 @@ OCL_TEST_P(ScaleAdd, Mat) //////////////////////////////// PatchNans //////////////////////////////////////////////// -PARAM_TEST_CASE(PatchNaNs, Channels, bool) +template +static _Tp randomNan(RNG& rng); + +template<> +float randomNan(RNG& rng) { + uint32_t r = rng.next(); + Cv32suf v; + v.u = r; + // set the exponent and one mantissa bit to get a NaN (and not an infinity) + v.u = v.u | 0x7f800001; + return v.f; +} + +template<> +double randomNan(RNG& rng) +{ + uint32_t r0 = rng.next(); + uint32_t r1 = rng.next(); + Cv64suf v; + v.u = (uint64_t(r0) << 32) | uint64_t(r1); + // set the exponent and one mantissa bit to get a NaN (and not an infinity) + v.u = v.u | 0x7ff0000000000001; + return v.f; +} + +PARAM_TEST_CASE(PatchNaNs, MatDepth, Channels, bool) +{ + int depth; int cn; bool use_roi; double value; @@ -1731,13 +1758,14 @@ PARAM_TEST_CASE(PatchNaNs, Channels, bool) virtual void SetUp() { - cn = GET_PARAM(0); - use_roi = GET_PARAM(1); + depth = GET_PARAM(0); + cn = GET_PARAM(1); + use_roi = GET_PARAM(2); } void generateTestData() { - const int type = CV_MAKE_TYPE(CV_32F, cn); + const int type = CV_MAKE_TYPE(depth, cn); Size roiSize = randomSize(1, 10); Border srcBorder = randomBorder(0, use_roi ? MAX_VALUE : 0); @@ -1747,9 +1775,15 @@ PARAM_TEST_CASE(PatchNaNs, Channels, bool) roiSize.width *= cn; for (int y = 0; y < roiSize.height; ++y) { - float * const ptr = src_roi.ptr(y); for (int x = 0; x < roiSize.width; ++x) - ptr[x] = randomInt(-1, 1) == 0 ? std::numeric_limits::quiet_NaN() : ptr[x]; + { + if (randomInt(-1, 1) != 0) + continue; + if (depth == CV_32F) + src_roi.ptr(y)[x] = randomNan(rng); + else + src_roi.ptr(y)[x] = randomNan(rng); + } } value = randomDouble(-100, 100); @@ -1950,7 +1984,7 @@ OCL_INSTANTIATE_TEST_CASE_P(Arithm, InRange, Combine(OCL_ALL_DEPTHS, OCL_ALL_CHA OCL_INSTANTIATE_TEST_CASE_P(Arithm, ConvertScaleAbs, Combine(OCL_ALL_DEPTHS, OCL_ALL_CHANNELS, Bool())); OCL_INSTANTIATE_TEST_CASE_P(Arithm, ConvertFp16, Combine(OCL_ALL_CHANNELS, Bool())); OCL_INSTANTIATE_TEST_CASE_P(Arithm, ScaleAdd, Combine(OCL_ALL_DEPTHS, OCL_ALL_CHANNELS, Bool())); -OCL_INSTANTIATE_TEST_CASE_P(Arithm, PatchNaNs, Combine(OCL_ALL_CHANNELS, Bool())); +OCL_INSTANTIATE_TEST_CASE_P(Arithm, PatchNaNs, Combine(testing::Values(CV_32F, CV_64F), OCL_ALL_CHANNELS, Bool())); OCL_INSTANTIATE_TEST_CASE_P(Arithm, Psnr, Combine(::testing::Values((MatDepth)CV_8U), OCL_ALL_CHANNELS, Bool())); OCL_INSTANTIATE_TEST_CASE_P(Arithm, UMatDot, Combine(OCL_ALL_DEPTHS, OCL_ALL_CHANNELS, Bool())); diff --git a/modules/core/test/test_math.cpp b/modules/core/test/test_math.cpp index 7b7f1b4a3f..8e5a47b6f7 100644 --- a/modules/core/test/test_math.cpp +++ b/modules/core/test/test_math.cpp @@ -4488,5 +4488,109 @@ testing::Values( INT_MIN) )); +/////////////////////////////////////////////////////////////////////////////////////////// + +typedef testing::TestWithParam Core_PatchNaNs; + +static float randomNanFlt(RNG& rng) +{ + uint32_t r = rng.next(); + Cv32suf v; + v.u = r; + // set the exponent and one mantissa bit to get a NaN (and not an infinity) + v.u = v.u | 0x7f800001; + return v.f; +} + +static double randomNanDbl(RNG& rng) +{ + uint32_t r0 = rng.next(); + uint32_t r1 = rng.next(); + Cv64suf v; + v.u = (uint64_t(r0) << 32) | uint64_t(r1); + // set the exponent and one mantissa bit to get a NaN (and not an infinity) + v.u = v.u | 0x7ff0000000000001; + return v.f; +} + +static void fillReferenceWithNans(cv::RNG& rng, const Size& sz, int type, double val, cv::Mat& in, cv::Mat& gold) +{ + in.create(sz, type); + gold.create(sz, type); + + for( int i = 0; i < in.rows; i++ ) + { + for( int j = 0; j < in.cols; j++ ) + { + if (CV_MAT_DEPTH(type) == CV_64F) + { + switch( rng.uniform(0, 4) ) + { + case 0: + in.at(i, j) = randomNanDbl(rng); + gold.at(i, j) = val; + break; + case 1: + in.at(i, j) = std::numeric_limits::infinity(); + gold.at(i, j) = in.at(i, j); + break; + case 2: + in.at(i, j) = -std::numeric_limits::infinity(); + gold.at(i, j) = in.at(i, j); + break; + + default: + in.at(i, j) = rng.uniform(-100.0, 100.0); + gold.at(i, j) = in.at(i, j); + break; + } + } + if (CV_MAT_DEPTH(type) == CV_32F) + { + switch( rng.uniform(0, 4) ) + { + case 0: + in.at(i, j) = randomNanFlt(rng); + gold.at(i, j) = val; + break; + case 1: + in.at(i, j) = std::numeric_limits::infinity(); + gold.at(i, j) = in.at(i, j); + break; + case 2: + in.at(i, j) = -std::numeric_limits::infinity(); + gold.at(i, j) = in.at(i, j); + break; + + default: + in.at(i, j) = rng.uniform(-100.0, 100.0); + gold.at(i, j) = in.at(i, j); + break; + } + } + } + } +} + +TEST_P(Core_PatchNaNs, accuracy) +{ + const int depth = GetParam(); + cv::Mat in, out, gold; + fillReferenceWithNans(theRNG(), cv::Size(127, 71), CV_MAKE_TYPE(depth, 1), 142., in, gold); + cv::patchNaNs(in, 142.); + // bit-exact check independant from + cv::Mat in_bin(in.rows, in.cols*in.elemSize(), CV_8UC1, in.data); + cv::Mat gold_bin(gold.rows, gold.cols*gold.elemSize(), CV_8UC1, gold.data); + EXPECT_EQ(0, cvtest::norm(gold_bin, in_bin, cv::NORM_INF)); +} + +TEST(Core_PatchNaNs_UnsupportedDepth, accuracy) +{ + Mat src(3, 3, CV_16F); + EXPECT_THROW(patchNaNs(src, 0), cv::Exception); +} + +INSTANTIATE_TEST_CASE_P(/* */, Core_PatchNaNs, testing::Values(CV_32F, CV_64F)); + }} // namespace /* End of file. */