From 0bbba847a413f86e91d81f28f9cc7bd04a3375c9 Mon Sep 17 00:00:00 2001 From: Andrey Kamaev Date: Sat, 15 Dec 2012 15:21:52 +0400 Subject: [PATCH 1/2] Fix equalization formula in equalizeHist function & rewrite in C++ Old implementation did lut[i] = 255 * (count(Y <= i)) / (width * height) which actually shifts uniform histograms. From now histogram is equalized as C = count(Y == min(Y)) lut[i] = 255 * (count(Y <= i) - C) / (width * height - C) --- modules/imgproc/src/histogram.cpp | 138 ++++++++++++++++++++---------- 1 file changed, 93 insertions(+), 45 deletions(-) diff --git a/modules/imgproc/src/histogram.cpp b/modules/imgproc/src/histogram.cpp index edcb24057..353ad5e0b 100644 --- a/modules/imgproc/src/histogram.cpp +++ b/modules/imgproc/src/histogram.cpp @@ -2407,58 +2407,106 @@ cvCalcProbDensity( const CvHistogram* hist, const CvHistogram* hist_mask, CV_IMPL void cvEqualizeHist( const CvArr* srcarr, CvArr* dstarr ) { - CvMat sstub, *src = cvGetMat(srcarr, &sstub); - CvMat dstub, *dst = cvGetMat(dstarr, &dstub); - - CV_Assert( CV_ARE_SIZES_EQ(src, dst) && CV_ARE_TYPES_EQ(src, dst) && - CV_MAT_TYPE(src->type) == CV_8UC1 ); - CvSize size = cvGetMatSize(src); - if( CV_IS_MAT_CONT(src->type & dst->type) ) - { - size.width *= size.height; - size.height = 1; - } - int x, y; - const int hist_sz = 256; - int hist[hist_sz]; - memset(hist, 0, sizeof(hist)); - - for( y = 0; y < size.height; y++ ) - { - const uchar* sptr = src->data.ptr + src->step*y; - for( x = 0; x < size.width; x++ ) - hist[sptr[x]]++; - } - - float scale = 255.f/(size.width*size.height); - int sum = 0; - uchar lut[hist_sz+1]; - - for( int i = 0; i < hist_sz; i++ ) - { - sum += hist[i]; - int val = cvRound(sum*scale); - lut[i] = CV_CAST_8U(val); - } - - lut[0] = 0; - for( y = 0; y < size.height; y++ ) - { - const uchar* sptr = src->data.ptr + src->step*y; - uchar* dptr = dst->data.ptr + dst->step*y; - for( x = 0; x < size.width; x++ ) - dptr[x] = lut[sptr[x]]; - } + cv::equalizeHist(cv::cvarrToMat(srcarr), cv::cvarrToMat(dstarr)); } - void cv::equalizeHist( InputArray _src, OutputArray _dst ) { Mat src = _src.getMat(); + CV_Assert( src.type() == CV_8UC1 ); + _dst.create( src.size(), src.type() ); Mat dst = _dst.getMat(); - CvMat _csrc = src, _cdst = dst; - cvEqualizeHist( &_csrc, &_cdst ); + + if(src.empty()) + return; + + const int hist_sz = (1 << (8*sizeof(uchar))); + int hist[hist_sz] = {0,}; + + const size_t sstep = src.step; + const size_t dstep = dst.step; + + int width = src.cols; + int height = src.rows; + + if (src.isContinuous()) + { + width *= height; + height = 1; + } + + for (const uchar* ptr = src.ptr(); height--; ptr += sstep) + { + int x = 0; + for (; x <= width - 4; x += 4) + { + int t0 = ptr[x], t1 = ptr[x+1]; + hist[t0]++; hist[t1]++; + t0 = ptr[x+2]; t1 = ptr[x+3]; + hist[t0]++; hist[t1]++; + } + + for (; x < width; ++x, ++ptr) + hist[ptr[x]]++; + } + + int i = 0; + while (!hist[i]) ++i; + + int total = (int)src.total(); + if (hist[i] == total) + { + dst.setTo(i); + return; + } + + float scale = (hist_sz - 1.f)/(total - hist[i]); + int sum = 0; + + int lut[hist_sz]; + + for (lut[i++] = 0; i < hist_sz; ++i) + { + sum += hist[i]; + lut[i] = saturate_cast(sum * scale); + } + + int cols = src.cols; + int rows = src.rows; + + if (src.isContinuous() && dst.isContinuous()) + { + cols *= rows; + rows = 1; + } + + const uchar* sptr = src.ptr(); + uchar* dptr = dst.ptr(); + + for (; rows--; sptr += sstep, dptr += dstep) + { + int x = 0; + for (; x <= cols - 4; x += 4) + { + int v0 = sptr[x]; + int v1 = sptr[x+1]; + int x0 = lut[v0]; + int x1 = lut[v1]; + dptr[x] = (uchar)x0; + dptr[x+1] = (uchar)x1; + + v0 = sptr[x+2]; + v1 = sptr[x+3]; + x0 = lut[v0]; + x1 = lut[v1]; + dptr[x+2] = (uchar)x0; + dptr[x+3] = (uchar)x1; + } + + for (; x < cols; ++x) + dptr[x] = (uchar)lut[sptr[x]]; + } } /* Implementation of RTTI and Generic Functions for CvHistogram */ From 98d7d992441c86055efea1cd508aa21ba3093346 Mon Sep 17 00:00:00 2001 From: Daniil Osokin Date: Tue, 18 Dec 2012 16:10:11 +0400 Subject: [PATCH 2/2] Add threaded version of equalizeHist --- modules/imgproc/src/histogram.cpp | 223 +++++++++++++++++++++--------- 1 file changed, 160 insertions(+), 63 deletions(-) diff --git a/modules/imgproc/src/histogram.cpp b/modules/imgproc/src/histogram.cpp index 353ad5e0b..12edfa0c2 100644 --- a/modules/imgproc/src/histogram.cpp +++ b/modules/imgproc/src/histogram.cpp @@ -2404,6 +2404,146 @@ cvCalcProbDensity( const CvHistogram* hist, const CvHistogram* hist_mask, } } +class EqualizeHistCalcHist_Invoker +{ +public: + enum {HIST_SZ = 256}; + +#ifdef HAVE_TBB + typedef tbb::mutex* MutextPtr; +#else + typedef void* MutextPtr; +#endif + + EqualizeHistCalcHist_Invoker(cv::Mat& src, int* histogram, MutextPtr histogramLock) + : src_(src), globalHistogram_(histogram), histogramLock_(histogramLock) + { } + + void operator()( const cv::BlockedRange& rowRange ) const + { + int localHistogram[HIST_SZ] = {0, }; + + const size_t sstep = src_.step; + + int width = src_.cols; + int height = rowRange.end() - rowRange.begin(); + + if (src_.isContinuous()) + { + width *= height; + height = 1; + } + + for (const uchar* ptr = src_.ptr(rowRange.begin()); height--; ptr += sstep) + { + int x = 0; + for (; x <= width - 4; x += 4) + { + int t0 = ptr[x], t1 = ptr[x+1]; + localHistogram[t0]++; localHistogram[t1]++; + t0 = ptr[x+2]; t1 = ptr[x+3]; + localHistogram[t0]++; localHistogram[t1]++; + } + + for (; x < width; ++x, ++ptr) + localHistogram[ptr[x]]++; + } + +#ifdef HAVE_TBB + tbb::mutex::scoped_lock lock(*histogramLock_); +#endif + + for( int i = 0; i < HIST_SZ; i++ ) + globalHistogram_[i] += localHistogram[i]; + } + + static bool isWorthParallel( const cv::Mat& src ) + { +#ifdef HAVE_TBB + return ( src.total() >= 640*480 ); +#else + (void)src; + return false; +#endif + } + +private: + EqualizeHistCalcHist_Invoker& operator=(const EqualizeHistCalcHist_Invoker&); + + cv::Mat& src_; + int* globalHistogram_; + MutextPtr histogramLock_; +}; + +class EqualizeHistLut_Invoker +{ +public: + EqualizeHistLut_Invoker( cv::Mat& src, cv::Mat& dst, int* lut ) + : src_(src), + dst_(dst), + lut_(lut) + { } + + void operator()( const cv::BlockedRange& rowRange ) const + { + const size_t sstep = src_.step; + const size_t dstep = dst_.step; + + int width = src_.cols; + int height = rowRange.end() - rowRange.begin(); + int* lut = lut_; + + if (src_.isContinuous() && dst_.isContinuous()) + { + width *= height; + height = 1; + } + + const uchar* sptr = src_.ptr(rowRange.begin()); + uchar* dptr = dst_.ptr(rowRange.begin()); + + for (; height--; sptr += sstep, dptr += dstep) + { + int x = 0; + for (; x <= width - 4; x += 4) + { + int v0 = sptr[x]; + int v1 = sptr[x+1]; + int x0 = lut[v0]; + int x1 = lut[v1]; + dptr[x] = (uchar)x0; + dptr[x+1] = (uchar)x1; + + v0 = sptr[x+2]; + v1 = sptr[x+3]; + x0 = lut[v0]; + x1 = lut[v1]; + dptr[x+2] = (uchar)x0; + dptr[x+3] = (uchar)x1; + } + + for (; x < width; ++x) + dptr[x] = (uchar)lut[sptr[x]]; + } + } + + static bool isWorthParallel( const cv::Mat& src ) + { +#ifdef HAVE_TBB + return ( src.total() >= 640*480 ); +#else + (void)src; + return false; +#endif + } + +private: + EqualizeHistLut_Invoker& operator=(const EqualizeHistLut_Invoker&); + + cv::Mat& src_; + cv::Mat& dst_; + int* lut_; +}; CV_IMPL void cvEqualizeHist( const CvArr* srcarr, CvArr* dstarr ) { @@ -2421,35 +2561,25 @@ void cv::equalizeHist( InputArray _src, OutputArray _dst ) if(src.empty()) return; - const int hist_sz = (1 << (8*sizeof(uchar))); +#ifdef HAVE_TBB + tbb::mutex histogramLockInstance; + EqualizeHistCalcHist_Invoker::MutextPtr histogramLock = &histogramLockInstance; +#else + EqualizeHistCalcHist_Invoker::MutextPtr histogramLock = 0; +#endif + + const int hist_sz = EqualizeHistCalcHist_Invoker::HIST_SZ; int hist[hist_sz] = {0,}; + int lut[hist_sz]; - const size_t sstep = src.step; - const size_t dstep = dst.step; + EqualizeHistCalcHist_Invoker calcBody(src, hist, histogramLock); + EqualizeHistLut_Invoker lutBody(src, dst, lut); + cv::BlockedRange heightRange(0, src.rows); - int width = src.cols; - int height = src.rows; - - if (src.isContinuous()) - { - width *= height; - height = 1; - } - - for (const uchar* ptr = src.ptr(); height--; ptr += sstep) - { - int x = 0; - for (; x <= width - 4; x += 4) - { - int t0 = ptr[x], t1 = ptr[x+1]; - hist[t0]++; hist[t1]++; - t0 = ptr[x+2]; t1 = ptr[x+3]; - hist[t0]++; hist[t1]++; - } - - for (; x < width; ++x, ++ptr) - hist[ptr[x]]++; - } + if(EqualizeHistCalcHist_Invoker::isWorthParallel(src)) + parallel_for(heightRange, calcBody); + else + calcBody(heightRange); int i = 0; while (!hist[i]) ++i; @@ -2464,49 +2594,16 @@ void cv::equalizeHist( InputArray _src, OutputArray _dst ) float scale = (hist_sz - 1.f)/(total - hist[i]); int sum = 0; - int lut[hist_sz]; - for (lut[i++] = 0; i < hist_sz; ++i) { sum += hist[i]; lut[i] = saturate_cast(sum * scale); } - int cols = src.cols; - int rows = src.rows; - - if (src.isContinuous() && dst.isContinuous()) - { - cols *= rows; - rows = 1; - } - - const uchar* sptr = src.ptr(); - uchar* dptr = dst.ptr(); - - for (; rows--; sptr += sstep, dptr += dstep) - { - int x = 0; - for (; x <= cols - 4; x += 4) - { - int v0 = sptr[x]; - int v1 = sptr[x+1]; - int x0 = lut[v0]; - int x1 = lut[v1]; - dptr[x] = (uchar)x0; - dptr[x+1] = (uchar)x1; - - v0 = sptr[x+2]; - v1 = sptr[x+3]; - x0 = lut[v0]; - x1 = lut[v1]; - dptr[x+2] = (uchar)x0; - dptr[x+3] = (uchar)x1; - } - - for (; x < cols; ++x) - dptr[x] = (uchar)lut[sptr[x]]; - } + if(EqualizeHistLut_Invoker::isWorthParallel(src)) + parallel_for(heightRange, lutBody); + else + lutBody(heightRange); } /* Implementation of RTTI and Generic Functions for CvHistogram */