Compare commits

...

2 Commits

Author SHA1 Message Date
Pratham Kumar
db0380a7a0 imgproc: parallelize calcHist using stripe-based parallel_for_ (#30585)
* imgproc: parallelize calcHist using stripe-based parallel_for_

* backport calcHist parallelization from 5.x
2026-11-09 12:21:40 +05:30
NIKHIL
54142329e9 imgproc: fix stackBlur results on full-range 16-bit images 🤖🤖🤖 (#30671)
* imgproc: fix stackBlur results on full-range 16-bit images

The 16-bit SIMD paths overflowed 32 bits in (sum * mul) >> shr and saturated when adding two 16-bit samples before the kernel multiply. Add the samples after widening, and divide as the scalar loops do, so SIMD and scalar agree and the result no longer depends on the thread count.

* imgproc: scale the stackBlur thread test with the machine's thread count

Restore the count with a guard so a throw can't leak it. Use an odd image width: with a width of 160, stripes at some thread counts, including 12, leave nothing to the scalar loop and the test passes on the unfixed code.

* imgproc: use 64-bit sums for large stackBlur kernels on 16-bit images

The largest sum is the largest absolute value times (radius + 1)^2, which no longer fits an int from radius 181 for CV_16U and radius 255 for CV_16S. Those radii now take the scalar loops with int64_t sums. Radii that fit keep the 32-bit SIMD paths.

* imgproc: use 64-bit sums for large stackBlur kernels on 8-bit images
2026-11-09 12:20:20 +05:30
4 changed files with 389 additions and 54 deletions

View File

@@ -806,6 +806,83 @@ namespace cv
}
#endif
namespace cv
{
template<typename CalcFunc>
static void calcHistRunParallel( std::vector<uchar*>& ptrs, const std::vector<int>& deltas,
Size imsize, Mat& hist, int dims, const float** ranges,
const double* uniranges, bool uniform, int esz1, CalcFunc calcFunc )
{
const int64 N = (int64)imsize.width * imsize.height;
const size_t M = hist.total();
if( N < (int64)256 * 1024 )
{
calcFunc( ptrs, deltas, imsize, hist, dims, ranges, uniranges, uniform );
return;
}
int nstripes = M > 0 ? std::max(1, std::min(8, cvRound((double)N * 1e-3 / (double)M))) : 1;
nstripes = std::min( nstripes, getNumThreads() );
if( nstripes <= 1 )
{
calcFunc( ptrs, deltas, imsize, hist, dims, ranges, uniranges, uniform );
return;
}
const bool splitByRows = imsize.height > 1;
const int splitDim = splitByRows ? imsize.height : imsize.width;
nstripes = std::min( nstripes, splitDim );
std::vector<Mat> local_hists( nstripes );
parallel_for_( Range( 0, nstripes ), [&]( const Range& r )
{
for( int s = r.start; s < r.end; s++ )
{
const int start = (int)((int64)s * splitDim / nstripes);
const int end = (int)((int64)(s + 1) * splitDim / nstripes);
if( start >= end )
continue;
std::vector<uchar*> local_ptrs = ptrs;
Size local_size;
if( splitByRows )
{
for( int i = 0; i < dims; i++ )
{
size_t row_step = ((size_t)imsize.width * deltas[i*2] + deltas[i*2 + 1]) * esz1;
local_ptrs[i] = (uchar*)((const uchar*)local_ptrs[i] + (size_t)start * row_step);
}
if( local_ptrs[dims] )
local_ptrs[dims] = (uchar*)((const uchar*)local_ptrs[dims] + (size_t)start * deltas[dims*2 + 1]);
local_size = Size( imsize.width, end - start );
}
else
{
for( int i = 0; i < dims; i++ )
local_ptrs[i] = (uchar*)((const uchar*)local_ptrs[i] + (size_t)start * deltas[i*2] * esz1);
if( local_ptrs[dims] )
local_ptrs[dims] = (uchar*)((const uchar*)local_ptrs[dims] + (size_t)start * deltas[dims*2]);
local_size = Size( end - start, 1 );
}
local_hists[s] = Mat::zeros( hist.dims, hist.size.p, hist.type() );
calcFunc( local_ptrs, deltas, local_size, local_hists[s], dims, ranges, uniranges, uniform );
}
});
for( int s = 0; s < nstripes; s++ )
{
if( !local_hists[s].empty() )
hist += local_hists[s];
}
}
} // namespace cv
void cv::calcHist( const Mat* images, int nimages, const int* channels,
InputArray _mask, OutputArray _hist, int dims, const int* histSize,
const float** ranges, bool uniform, bool accumulate )
@@ -857,13 +934,14 @@ void cv::calcHist( const Mat* images, int nimages, const int* channels,
const double* _uniranges = uniform ? &uniranges[0] : 0;
int depth = images[0].depth();
int esz1 = (int)images[0].elemSize1();
if( depth == CV_8U )
calcHist_8u(ptrs, deltas, imsize, ihist, dims, ranges, _uniranges, uniform );
calcHistRunParallel( ptrs, deltas, imsize, ihist, dims, ranges, _uniranges, uniform, esz1, &calcHist_8u );
else if( depth == CV_16U )
calcHist_<ushort>(ptrs, deltas, imsize, ihist, dims, ranges, _uniranges, uniform );
calcHistRunParallel( ptrs, deltas, imsize, ihist, dims, ranges, _uniranges, uniform, esz1, &calcHist_<ushort> );
else if( depth == CV_32F )
calcHist_<float>(ptrs, deltas, imsize, ihist, dims, ranges, _uniranges, uniform );
calcHistRunParallel( ptrs, deltas, imsize, ihist, dims, ranges, _uniranges, uniform, esz1, &calcHist_<float> );
else
CV_Error(cv::Error::StsUnsupportedFormat, "");

View File

@@ -42,6 +42,8 @@ OTHER DEALINGS IN THE SOFTWARE.
#include <iostream>
#include <algorithm>
#include <climits>
#include <limits>
using namespace std;
@@ -96,6 +98,12 @@ inline int opRow(const T* , T* , const std::vector<ushort>& , const float , cons
return radius * CN;
}
// A 16-bit sum times stackblurMul doesn't fit in 32 bits, so 16-bit data divides as the scalar loops do.
inline v_int32 stackblurDivide(const v_int32& sum, const v_float32& mulVal)
{
return v_trunc(v_mul(v_cvt_f32(sum), mulVal));
}
template<>
inline int opRow<uchar>(const uchar* srcPtr, uchar* dstPtr, const std::vector<ushort>& kVec, const float , const int radius, const int CN, const int widthCN)
{
@@ -195,7 +203,7 @@ inline int opRow<uchar>(const uchar* srcPtr, uchar* dstPtr, const std::vector<us
}
template<>
inline int opRow<ushort>(const ushort* srcPtr, ushort* dstPtr, const std::vector<ushort>& kVec, const float , const int radius, const int CN, const int widthCN)
inline int opRow<ushort>(const ushort* srcPtr, ushort* dstPtr, const std::vector<ushort>& kVec, const float mulVal, const int radius, const int CN, const int widthCN)
{
int kernelSize = (int)kVec.size();
@@ -203,12 +211,9 @@ inline int opRow<ushort>(const ushort* srcPtr, ushort* dstPtr, const std::vector
if (radius > STACKBLUR_MAX_RADIUS)
return i;
const int mulValTab= stackblurMul[radius];
const int shrValTab= stackblurShr[radius];
const int VEC_LINE = VTraits<v_uint16>::vlanes();
v_uint32 v_mulVal = vx_setall_u32(mulValTab);
v_float32 v_mulVal = vx_setall_f32(mulVal);
if (kernelSize == 3)
{
for (; i <= widthCN - VEC_LINE; i += VEC_LINE)
@@ -221,7 +226,7 @@ inline int opRow<ushort>(const ushort* srcPtr, ushort* dstPtr, const std::vector
x1l = v_add(v_add(x1l, x1l), v_add(x0l, x2l));
x1h = v_add(v_add(x1h, x1h), v_add(x0h, x2h));
v_store(dstPtr + i, v_pack(v_shr(v_mul(x1l, v_mulVal), shrValTab), v_shr(v_mul(x1h, v_mulVal), shrValTab)));
v_store(dstPtr + i, v_pack_u(stackblurDivide(v_reinterpret_as_s32(x1l), v_mulVal), stackblurDivide(v_reinterpret_as_s32(x1h), v_mulVal)));
}
}
else
@@ -243,11 +248,18 @@ inline int opRow<ushort>(const ushort* srcPtr, ushort* dstPtr, const std::vector
v_uint16 k1 = vx_setall_u16(kx[k]);
v_uint16 k2 = vx_setall_u16(kx[k + 1]);
// Adding two 16-bit values first would saturate, so each one is multiplied on its own.
v_uint32 y0, y1;
v_mul_expand(v_add(vx_load(srcPtr - j), vx_load(srcPtr + j)), k1, y0, y1);
v_mul_expand(vx_load(srcPtr - j), k1, y0, y1);
s0 = v_add(s0, y0);
s1 = v_add(s1, y1);
v_mul_expand(v_add(vx_load(srcPtr - j - CN), vx_load(srcPtr + j + CN)), k2, y0, y1);
v_mul_expand(vx_load(srcPtr + j), k1, y0, y1);
s0 = v_add(s0, y0);
s1 = v_add(s1, y1);
v_mul_expand(vx_load(srcPtr - j - CN), k2, y0, y1);
s0 = v_add(s0, y0);
s1 = v_add(s1, y1);
v_mul_expand(vx_load(srcPtr + j + CN), k2, y0, y1);
s0 = v_add(s0, y0);
s1 = v_add(s1, y1);
}
@@ -256,15 +268,15 @@ inline int opRow<ushort>(const ushort* srcPtr, ushort* dstPtr, const std::vector
v_uint16 k1 = vx_setall_u16(kx[k]);
v_uint32 y0, y1;
v_mul_expand(v_add(vx_load(srcPtr - j), vx_load(srcPtr + j)), k1, y0, y1);
v_mul_expand(vx_load(srcPtr - j), k1, y0, y1);
s0 = v_add(s0, y0);
s1 = v_add(s1, y1);
v_mul_expand(vx_load(srcPtr + j), k1, y0, y1);
s0 = v_add(s0, y0);
s1 = v_add(s1, y1);
}
s0 = v_shr(v_mul(s0, v_mulVal), shrValTab);
s1 = v_shr(v_mul(s1, v_mulVal), shrValTab);
v_store(dstPtr + i, v_pack(s0, s1));
v_store(dstPtr + i, v_pack_u(stackblurDivide(v_reinterpret_as_s32(s0), v_mulVal), stackblurDivide(v_reinterpret_as_s32(s1), v_mulVal)));
}
}
@@ -272,7 +284,7 @@ inline int opRow<ushort>(const ushort* srcPtr, ushort* dstPtr, const std::vector
}
template<>
inline int opRow<short>(const short* srcPtr, short* dstPtr, const std::vector<ushort>& kVec, const float , const int radius, const int CN, const int widthCN)
inline int opRow<short>(const short* srcPtr, short* dstPtr, const std::vector<ushort>& kVec, const float mulVal, const int radius, const int CN, const int widthCN)
{
int kernelSize = (int)kVec.size();
int i = radius * CN;
@@ -280,11 +292,8 @@ inline int opRow<short>(const short* srcPtr, short* dstPtr, const std::vector<us
if (radius > STACKBLUR_MAX_RADIUS)
return i;
const int mulValTab= stackblurMul[radius];
const int shrValTab= stackblurShr[radius];
const int VEC_LINE = VTraits<v_int16>::vlanes();
v_int32 v_mulVal = vx_setall_s32(mulValTab);
v_float32 v_mulVal = vx_setall_f32(mulVal);
if (kernelSize == 3)
{
@@ -298,7 +307,7 @@ inline int opRow<short>(const short* srcPtr, short* dstPtr, const std::vector<us
x1l = v_add(v_add(x1l, x1l), v_add(x0l, x2l));
x1h = v_add(v_add(x1h, x1h), v_add(x0h, x2h));
v_store(dstPtr + i, v_pack(v_shr(v_mul(x1l, v_mulVal), shrValTab), v_shr(v_mul(x1h, v_mulVal), shrValTab)));
v_store(dstPtr + i, v_pack(stackblurDivide(x1l, v_mulVal), stackblurDivide(x1h, v_mulVal)));
}
}
else
@@ -319,12 +328,19 @@ inline int opRow<short>(const short* srcPtr, short* dstPtr, const std::vector<us
v_int16 k1 = vx_setall_s16((short)kx[k]);
v_int16 k2 = vx_setall_s16((short)kx[k + 1]);
// Adding two 16-bit values first would saturate, so each one is multiplied on its own.
v_int32 y0, y1;
v_mul_expand(v_add(vx_load(srcPtr - j), vx_load(srcPtr + j)), k1, y0, y1);
v_mul_expand(vx_load(srcPtr - j), k1, y0, y1);
s0 = v_add(s0, y0);
s1 = v_add(s1, y1);
v_mul_expand(v_add(vx_load(srcPtr - j - CN), vx_load(srcPtr + j + CN)), k2, y0, y1);
v_mul_expand(vx_load(srcPtr + j), k1, y0, y1);
s0 = v_add(s0, y0);
s1 = v_add(s1, y1);
v_mul_expand(vx_load(srcPtr - j - CN), k2, y0, y1);
s0 = v_add(s0, y0);
s1 = v_add(s1, y1);
v_mul_expand(vx_load(srcPtr + j + CN), k2, y0, y1);
s0 = v_add(s0, y0);
s1 = v_add(s1, y1);
}
@@ -332,15 +348,15 @@ inline int opRow<short>(const short* srcPtr, short* dstPtr, const std::vector<us
{
v_int16 k1 = vx_setall_s16((short)kx[k]);
v_int32 y0, y1;
v_mul_expand(v_add(vx_load(srcPtr - j), vx_load(srcPtr + j)), k1, y0, y1);
v_mul_expand(vx_load(srcPtr - j), k1, y0, y1);
s0 = v_add(s0, y0);
s1 = v_add(s1, y1);
v_mul_expand(vx_load(srcPtr + j), k1, y0, y1);
s0 = v_add(s0, y0);
s1 = v_add(s1, y1);
}
s0 = v_shr(v_mul(s0, v_mulVal), shrValTab);
s1 = v_shr(v_mul(s1, v_mulVal), shrValTab);
v_store(dstPtr + i, v_pack(s0, s1));
v_store(dstPtr + i, v_pack(stackblurDivide(s0, v_mulVal), stackblurDivide(s1, v_mulVal)));
}
}
return i;
@@ -456,7 +472,7 @@ public:
{
width= dst.cols;
wm = width - 1;
mulVal = 1.0f / ((radius + 1) * (radius + 1));
mulVal = 1.0f / ((float)(radius + 1) * (radius + 1));
CN = src.channels();
}
@@ -565,9 +581,9 @@ public:
const T* srcPtr = src.ptr<T>(row);
T* dstPtr = dst.ptr<T>(row);
int radiusMul = (radius + 2) * (radius + 1) / 2;
const int64_t radiusMul = (int64_t)(radius + 2) * (radius + 1) / 2;
for (int ci = 0; ci < CN; ci++)
sum[ci] += (TBuf)srcPtr[ci] * radiusMul;
sum[ci] += (TBuf)srcPtr[ci] * (TBuf)radiusMul;
// compute diff
const T* srcPtr0 = srcPtr;
@@ -580,7 +596,7 @@ public:
{
diff[i*CN + ci] = (TBuf)srcPtr0[ci] - (TBuf)srcPtr[ci];
diffVal[ci] += diff[i*CN + ci];
sum[ci] += srcPtr0[ci] * (radius - i);
sum[ci] += (TBuf)srcPtr0[ci] * (radius - i);
}
}
@@ -917,7 +933,7 @@ inline int opColumn<uchar, int>(const uchar* srcPtr, uchar* dstPtr, uchar* stack
template<>
inline int opColumn<short, int>(const short* srcPtr, short* dstPtr, short* stack, int* sum, int* sumIn,
int* sumOut, const float , const int mulValTab, const int shrValTab,
int* sumOut, const float mulVal, const int mulValTab, const int shrValTab,
const int widthLen, const int ss, const int sp1)
{
int k = 0;
@@ -925,7 +941,7 @@ inline int opColumn<short, int>(const short* srcPtr, short* dstPtr, short* stack
{
const int VEC_LINE_16 = VTraits<v_int16>::vlanes();
const int VEC_LINE_32 = VTraits<v_int32>::vlanes();
v_int32 v_mulVal = vx_setall_s32(mulValTab);
v_float32 v_mulVal = vx_setall_f32(mulVal);
auto stackStartPtr = stack + ss * widthLen;
auto stackSp1Ptr = stack + sp1 * widthLen;
@@ -944,7 +960,7 @@ inline int opColumn<short, int>(const short* srcPtr, short* dstPtr, short* stack
v_sumOut0 = vx_load(sumOut + k);
v_sumOut1 = vx_load(sumOut + k + VEC_LINE_32);
v_store(dstPtr + k,v_pack(v_shr(v_mul(v_sum0, v_mulVal), shrValTab), v_shr(v_mul(v_sum1, v_mulVal), shrValTab)));
v_store(dstPtr + k, v_pack(stackblurDivide(v_sum0, v_mulVal), stackblurDivide(v_sum1, v_mulVal)));
v_sum0 = v_sub(v_sum0, v_sumOut0);
v_sum1 = v_sub(v_sum1, v_sumOut1);
@@ -987,7 +1003,7 @@ inline int opColumn<short, int>(const short* srcPtr, short* dstPtr, short* stack
template<>
inline int opColumn<ushort, int>(const ushort* srcPtr, ushort* dstPtr, ushort* stack, int* sum, int* sumIn,
int* sumOut, const float , const int mulValTab, const int shrValTab,
int* sumOut, const float mulVal, const int mulValTab, const int shrValTab,
const int widthLen, const int ss, const int sp1)
{
int k = 0;
@@ -995,7 +1011,7 @@ inline int opColumn<ushort, int>(const ushort* srcPtr, ushort* dstPtr, ushort* s
{
const int VEC_LINE_16 = VTraits<v_uint16>::vlanes();
const int VEC_LINE_32 = VTraits<v_int32>::vlanes();
v_uint32 v_mulVal = vx_setall_u32((uint32_t)mulValTab);
v_float32 v_mulVal = vx_setall_f32(mulVal);
auto stackStartPtr = stack + ss * widthLen;
auto stackSp1Ptr = stack + sp1 * widthLen;
@@ -1014,7 +1030,7 @@ inline int opColumn<ushort, int>(const ushort* srcPtr, ushort* dstPtr, ushort* s
v_sumOut0 = vx_load(sumOut + k);
v_sumOut1 = vx_load(sumOut + k + VEC_LINE_32);
v_store(dstPtr + k, v_pack(v_shr(v_mul(v_reinterpret_as_u32(v_sum0), v_mulVal), shrValTab), v_shr(v_mul(v_reinterpret_as_u32(v_sum1), v_mulVal), shrValTab)));
v_store(dstPtr + k, v_pack_u(stackblurDivide(v_sum0, v_mulVal), stackblurDivide(v_sum1, v_mulVal)));
v_sum0 = v_sub(v_sum0, v_sumOut0);
v_sum1 = v_sub(v_sum1, v_sumOut1);
@@ -1068,7 +1084,7 @@ public:
widthElem = CN * src.cols;
height = src.rows;
hm = src.rows - 1;
mulVal = 1.0f / ((radius + 1)*(radius + 1));
mulVal = 1.0f / ((float)(radius + 1) * (radius + 1));
if (radius <= STACKBLUR_MAX_RADIUS)
{
shrValTab = stackblurShr[radius];
@@ -1111,7 +1127,7 @@ public:
for (int k = 0; k < widthLen; k++)
{
stack[i * widthLen + k] = *(srcPtr + k);
sum[k] += *(srcPtr + k) * (i + 1);
sum[k] += (TBuf)*(srcPtr + k) * (i + 1);
sumOut[k] += *(srcPtr + k);
}
}
@@ -1123,7 +1139,7 @@ public:
{
T tmp = *(srcPtr + k);
stack[(i + radius) * widthLen + k] = tmp;
sum[k] += tmp * (radius - i + 1);
sum[k] += (TBuf)tmp * (radius - i + 1);
sumIn[k] += tmp;
}
}
@@ -1192,6 +1208,34 @@ private:
int shrValTab;
};
// The largest sum is the largest absolute value times (radius + 1)^2.
template<typename T>
static inline bool stackBlurSumFitsInt(int radius)
{
const int64_t maxAbs = std::max(-(int64_t)std::numeric_limits<T>::min(), (int64_t)std::numeric_limits<T>::max());
return (int64_t)(radius + 1) * (radius + 1) <= INT_MAX / maxAbs;
}
// The SIMD paths need 32-bit sums, so larger kernels take the scalar loops with 64-bit sums.
template<typename T>
static void stackBlurInteger(const Mat& src, Mat& dst, int radiusW, int radiusH, int widthElem, int numOfThreads)
{
if (radiusW != 0)
{
if (stackBlurSumFitsInt<T>(radiusW))
parallel_for_(Range(0, src.rows), ParallelStackBlurRow<T, int>(src, dst, radiusW), numOfThreads);
else
parallel_for_(Range(0, src.rows), ParallelStackBlurRow<T, int64_t>(src, dst, radiusW), numOfThreads);
}
if (radiusH != 0)
{
if (stackBlurSumFitsInt<T>(radiusH))
parallel_for_(Range(0, widthElem), ParallelStackBlurColumn<T, int>(dst, dst, radiusH), numOfThreads);
else
parallel_for_(Range(0, widthElem), ParallelStackBlurColumn<T, int64_t>(dst, dst, radiusH), numOfThreads);
}
}
void stackBlur(InputArray _src, OutputArray _dst, Size ksize)
{
CV_INSTRUMENT_REGION();
@@ -1232,24 +1276,15 @@ void stackBlur(InputArray _src, OutputArray _dst, Size ksize)
if (sdepth == CV_8U)
{
if (ksize.width != 1)
parallel_for_(Range(0, src.rows), ParallelStackBlurRow<uchar, int>(src, dst, radiusW), numOfThreads);
if (ksize.height != 1)
parallel_for_(Range(0, widthElem), ParallelStackBlurColumn<uchar, int>(dst, dst, radiusH), numOfThreads);
stackBlurInteger<uchar>(src, dst, radiusW, radiusH, widthElem, numOfThreads);
}
else if (sdepth == CV_16S)
{
if (ksize.width != 1)
parallel_for_(Range(0, src.rows), ParallelStackBlurRow<short, int>(src, dst, radiusW), numOfThreads);
if (ksize.height != 1)
parallel_for_(Range(0, widthElem), ParallelStackBlurColumn<short, int>(dst, dst, radiusH), numOfThreads);
stackBlurInteger<short>(src, dst, radiusW, radiusH, widthElem, numOfThreads);
}
else if (sdepth == CV_16U)
{
if (ksize.width != 1)
parallel_for_(Range(0, src.rows), ParallelStackBlurRow<ushort, int>(src, dst, radiusW), numOfThreads);
if (ksize.height != 1)
parallel_for_(Range(0, widthElem), ParallelStackBlurColumn<ushort, int>(dst, dst, radiusH), numOfThreads);
stackBlurInteger<ushort>(src, dst, radiusW, radiusH, widthElem, numOfThreads);
}
else if (sdepth == CV_32F)
{

View File

@@ -2032,6 +2032,73 @@ TEST(Imgproc_Hist_Calc, IPP_ranges_with_nonequal_exponent_21595)
ASSERT_EQ(histogram_u.at<float>(2), 4.f) << "1 not counts correctly, res: " << histogram_u.at<float>(2);
}
typedef testing::TestWithParam< testing::tuple<int, int, Size> > Imgproc_Hist_Calc_Parallel;
TEST_P(Imgproc_Hist_Calc_Parallel, consistency)
{
int type = testing::get<0>(GetParam());
int dims = testing::get<1>(GetParam());
Size sz = testing::get<2>(GetParam());
int cn = dims;
Mat buf(Size(sz.width + 1, sz.height), CV_MAKETYPE(type, cn));
randu(buf, Scalar::all(0), Scalar::all(256));
Mat src = buf(Rect(0, 0, sz.width, sz.height));
Mat mask;
if( dims == 1 )
{
mask.create(sz, CV_8UC1);
randu(mask, Scalar::all(0), Scalar::all(2));
}
std::vector<int> ch(dims), hs(dims);
std::vector<std::vector<float> > rv(dims);
std::vector<const float*> rp(dims);
for( int i = 0; i < dims; i++ )
{
ch[i] = i;
hs[i] = (dims == 1) ? 64 : (dims == 2) ? 8 : 4;
rv[i].assign(2, 0.f);
rv[i][1] = 256.f;
rp[i] = &rv[i][0];
}
// Row-by-row: every call stays below the parallel threshold.
Mat baseline;
for( int y = 0; y < src.rows; y++ )
{
Mat row = src.row(y);
Mat mask_row = (dims == 1) ? mask.row(y) : Mat();
bool acc = (y > 0);
calcHist(&row, 1, ch.data(), mask_row, baseline, dims, hs.data(), rp.data(), true, acc);
}
setNumThreads(8);
Mat parallel_result;
calcHist(&src, 1, ch.data(), mask, parallel_result, dims, hs.data(), rp.data(), true, false);
// Flattened to one row: column-split path.
Mat contiguous = src.clone();
Mat contiguous_result;
calcHist(&contiguous, 1, ch.data(), mask, contiguous_result, dims, hs.data(), rp.data(), true, false);
setNumThreads(-1);
EXPECT_EQ(0.0, cv::norm(baseline, parallel_result, NORM_INF))
<< "row-split, type=" << type << " dims=" << dims << " size=" << sz;
EXPECT_EQ(0.0, cv::norm(baseline, contiguous_result, NORM_INF))
<< "column-split, type=" << type << " dims=" << dims << " size=" << sz;
}
INSTANTIATE_TEST_CASE_P(/*nothing*/, Imgproc_Hist_Calc_Parallel,
testing::Combine(
testing::Values(CV_8U, CV_16U, CV_32F),
testing::Values(1, 2, 3),
testing::Values(Size(640, 480), Size(1920, 1080), Size(7, 13))
)
);
////////////////////////////////////////// equalizeHist() /////////////////////////////////////////
void equalizeHistReference(const Mat& src, Mat& dst)

View File

@@ -324,5 +324,160 @@ TEST(Imgproc_StackBlur, regression_28233)
EXPECT_NO_THROW(stackBlur(src2, dst2, Size(11, 11)));
EXPECT_EQ(dst2.at<uchar>(1, 1), 50);
}
typedef testing::TestWithParam<tuple<int, int, int> > StackBlur_FullRange;
// A constant image stays constant, however large its values, and full-range noise matches the reference.
TEST_P(StackBlur_FullRange, regression)
{
const int depth = get<0>(GetParam());
const int cn = get<1>(GetParam());
const int kernel = get<2>(GetParam());
const Size ksize(kernel, kernel);
const double maxVal = depth == CV_16U ? 65535 : 32767;
const double minVal = depth == CV_16U ? 0 : -32768;
for (double value : {maxVal, minVal, maxVal / 2})
{
Mat constant(64, 96, CV_MAKETYPE(depth, cn), Scalar::all(value)), dst;
stackBlur(constant, dst, ksize);
EXPECT_LE(cvtest::norm(constant, dst, NORM_INF), 1.) << "value=" << value;
}
// stackBlurRef accumulates in float, which drifts at this range, so the reference is the
// triangular kernel stackBlur stands for, computed in double.
const int radius = kernel / 2;
Mat tri(kernel, 1, CV_64F);
for (int i = 0; i < kernel; i++)
tri.at<double>(i) = (radius + 1 - std::abs(i - radius)) / double((radius + 1) * (radius + 1));
Mat noise(64, 96, CV_MAKETYPE(depth, cn)), dst, ref, dst64;
randu(noise, Scalar::all(minVal), Scalar::all(maxVal));
sepFilter2D(noise, ref, CV_64F, tri, tri, Point(-1, -1), 0, BORDER_REPLICATE);
stackBlur(noise, dst, ksize);
dst.convertTo(dst64, CV_64F);
const Rect inner(radius, radius, noise.cols - 2 * radius, noise.rows - 2 * radius);
EXPECT_LE(cvtest::norm(ref(inner), dst64(inner), NORM_INF), 2.);
}
INSTANTIATE_TEST_CASE_P(Imgproc, StackBlur_FullRange,
testing::Combine(
testing::Values(CV_16U, CV_16S),
testing::Values(1, 3, 4),
testing::Values(3, 5, 7, 15, 51)
)
);
typedef testing::TestWithParam<tuple<int, int> > StackBlur_LargeKernel;
// The sum no longer fits 32 bits from radius 181 for CV_16U and radius 255 for CV_16S.
TEST_P(StackBlur_LargeKernel, regression)
{
const int depth = get<0>(GetParam());
const int kernel = get<1>(GetParam());
const double maxVal = depth == CV_16U ? 65535 : 32767;
const double minVal = depth == CV_16U ? 0 : -32768;
const int side = kernel + 32;
for (double value : {maxVal, minVal, maxVal / 2})
{
Mat constant(side, side, CV_MAKETYPE(depth, 3), Scalar::all(value)), dst;
for (const Size& ksize : {Size(kernel, kernel), Size(kernel, 1), Size(1, kernel)})
{
stackBlur(constant, dst, ksize);
EXPECT_LE(cvtest::norm(constant, dst, NORM_INF), 1.) << "value=" << value << " ksize=" << ksize;
}
}
const int radius = kernel / 2;
Mat tri(kernel, 1, CV_64F);
for (int i = 0; i < kernel; i++)
tri.at<double>(i) = (radius + 1 - std::abs(i - radius)) / double((radius + 1) * (radius + 1));
Mat noise(side, side, depth), dst, ref, dst64;
randu(noise, Scalar::all(minVal), Scalar::all(maxVal));
sepFilter2D(noise, ref, CV_64F, tri, tri, Point(-1, -1), 0, BORDER_REPLICATE);
stackBlur(noise, dst, Size(kernel, kernel));
dst.convertTo(dst64, CV_64F);
const Rect inner(radius, radius, noise.cols - 2 * radius, noise.rows - 2 * radius);
EXPECT_LE(cvtest::norm(ref(inner), dst64(inner), NORM_INF), 2.);
}
INSTANTIATE_TEST_CASE_P(Imgproc, StackBlur_LargeKernel,
testing::Combine(
testing::Values(CV_16U, CV_16S),
testing::Values(361, 363, 509, 511, 513)
)
);
// For CV_8U the sum no longer fits 32 bits from radius 2901.
TEST(Imgproc_StackBlur, large_kernel_8u)
{
for (int kernel : {5801, 5803})
{
const int radius = kernel / 2;
for (double value : {255., 128.})
{
Mat row(3, kernel + 32, CV_8UC1, Scalar::all(value)), col(kernel + 32, 3, CV_8UC1, Scalar::all(value)), dst;
stackBlur(row, dst, Size(kernel, 1));
EXPECT_LE(cvtest::norm(row, dst, NORM_INF), 1.) << "value=" << value << " kernel=" << kernel;
stackBlur(col, dst, Size(1, kernel));
EXPECT_LE(cvtest::norm(col, dst, NORM_INF), 1.) << "value=" << value << " kernel=" << kernel;
}
Mat tri(1, kernel, CV_64F), one(1, 1, CV_64F, Scalar::all(1.));
for (int i = 0; i < kernel; i++)
tri.at<double>(i) = (radius + 1 - std::abs(i - radius)) / double((radius + 1) * (radius + 1));
Mat noise(3, 2 * kernel, CV_8UC1), dst, ref, dst64;
randu(noise, Scalar::all(0), Scalar::all(256));
sepFilter2D(noise, ref, CV_64F, tri, one, Point(-1, -1), 0, BORDER_REPLICATE);
stackBlur(noise, dst, Size(kernel, 1));
dst.convertTo(dst64, CV_64F);
const Rect inner(radius, 0, noise.cols - 2 * radius, noise.rows);
EXPECT_LE(cvtest::norm(ref(inner), dst64(inner), NORM_INF), 1.) << "kernel=" << kernel;
}
}
// (radius + 1)^2 itself doesn't fit 32 bits here.
TEST(Imgproc_StackBlur, full_range_huge_kernel)
{
const int kernel = 92681;
for (int depth : {CV_16U, CV_16S})
{
const double value = depth == CV_16U ? 65535 : -32768;
Mat constant(1, kernel + 2, depth, Scalar::all(value)), dst;
stackBlur(constant, dst, Size(kernel, 1));
EXPECT_LE(cvtest::norm(constant, dst, NORM_INF), 1.) << "depth=" << depth;
}
}
class ScopedThreadsGuard
{
public:
ScopedThreadsGuard() : old_threads(getNumThreads()) {}
~ScopedThreadsGuard() { setNumThreads(old_threads); }
void set(int n) { setNumThreads(n); }
private:
int old_threads;
};
// The result doesn't depend on how the image is split between threads. The width is
// odd so that some stripe leaves pixels to the scalar loop at any thread count.
TEST(Imgproc_StackBlur, full_range_thread_invariance)
{
Mat noise(120, 161, CV_16UC3), single, multi;
randu(noise, Scalar::all(0), Scalar::all(65535));
ScopedThreadsGuard threadsGuard;
const int threads = std::max(2, getNumThreads());
threadsGuard.set(1);
stackBlur(noise, single, Size(5, 5));
threadsGuard.set(threads);
stackBlur(noise, multi, Size(5, 5));
EXPECT_EQ(cvtest::norm(single, multi, NORM_INF), 0.);
}
}
}