Merge pull request #226 from taka-no-me/equalizeHist

Fix equalization formula in equalizeHist function & rewrite in C++

modules/imgproc/src/histogram.cpp

@@ -2404,61 +2404,206 @@ cvCalcProbDensity( const CvHistogram* hist, const CvHistogram* hist_mask,
     }
 }
 
-
-CV_IMPL void cvEqualizeHist( const CvArr* srcarr, CvArr* dstarr )
+class EqualizeHistCalcHist_Invoker
 {
-    CvMat sstub, *src = cvGetMat(srcarr, &sstub);
-    CvMat dstub, *dst = cvGetMat(dstarr, &dstub);
+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)
+    { }
 
-    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) )
+    void operator()( const cv::BlockedRange& rowRange ) const
     {
-        size.width *= size.height;
-        size.height = 1;
+        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<uchar>(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];
     }
-    int x, y;
-    const int hist_sz = 256;
-    int hist[hist_sz];
-    memset(hist, 0, sizeof(hist));
 
-    for( y = 0; y < size.height; y++ )
+    static bool isWorthParallel( const cv::Mat& src )
     {
-        const uchar* sptr = src->data.ptr + src->step*y;
-        for( x = 0; x < size.width; x++ )
-            hist[sptr[x]]++;
+#ifdef HAVE_TBB
+        return ( src.total() >= 640*480 );
+#else
+        (void)src;
+        return false;
+#endif
     }
 
-    float scale = 255.f/(size.width*size.height);
-    int sum = 0;
-    uchar lut[hist_sz+1];
+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)
+    { }
 
-    for( int i = 0; i < hist_sz; i++ )
+    void operator()( const cv::BlockedRange& rowRange ) const
     {
-        sum += hist[i];
-        int val = cvRound(sum*scale);
-        lut[i] = CV_CAST_8U(val);
+        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<uchar>(rowRange.begin());
+        uchar* dptr = dst_.ptr<uchar>(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]];
+        }
     }
 
-    lut[0] = 0;
-    for( y = 0; y < size.height; y++ )
+    static bool isWorthParallel( const cv::Mat& src )
     {
-        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]];
+#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 )
+{
+    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;
+
+#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];
+
+    EqualizeHistCalcHist_Invoker calcBody(src, hist, histogramLock);
+    EqualizeHistLut_Invoker      lutBody(src, dst, lut);
+    cv::BlockedRange heightRange(0, src.rows);
+
+    if(EqualizeHistCalcHist_Invoker::isWorthParallel(src))
+        parallel_for(heightRange, calcBody);
+    else
+        calcBody(heightRange);
+
+    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;
+
+    for (lut[i++] = 0; i < hist_sz; ++i)
+    {
+        sum += hist[i];
+        lut[i] = saturate_cast<uchar>(sum * scale);
+    }
+
+    if(EqualizeHistLut_Invoker::isWorthParallel(src))
+        parallel_for(heightRange, lutBody);
+    else
+        lutBody(heightRange);
 }
 
 /* Implementation of RTTI and Generic Functions for CvHistogram */