OpenCV源码解析：直方图均衡化的详细算法和过程

柔光的暖阳◎ 2022-05-16 03:42 198阅读 0赞

## 理论基础 ##

所谓直方图，在图像中，指的就是各个像素的统计值，就是一个像素在整幅图像中出现次数。

例如下面这张16个像素的图片，其直方图就是

![70][]![70 1][]

直方图均衡化，是将给定图像的直方图改造成均匀分布的直方图，从而扩大像素灰度值的动态范围，达到增强图像对比度的效果。

![70 2][]![70 3][]

![70 4][]![70 5][]

## OpenCV中的直方图均衡化 ##

OpneCv中，可以用calcHist进行图像的均衡化，也可以使用equalizeHist可以进行全局直方图均衡化（就是直接对整个图像的直方图进行均衡化）。这里以equalizeHist为例进行详细讲解。

## 先看一下equalizeHist的源码 ##

void cv::equalizeHist( InputArray _src, OutputArray _dst )
    {
        CV_INSTRUMENT_REGION()
    
        CV_Assert( _src.type() == CV_8UC1 );
    
        if (_src.empty())
            return;
    
        CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
                   ocl_equalizeHist(_src, _dst))
    
        Mat src = _src.getMat();
        _dst.create( src.size(), src.type() );
        Mat dst = _dst.getMat();
    
        CV_OVX_RUN(!ovx::skipSmallImages<VX_KERNEL_EQUALIZE_HISTOGRAM>(src.cols, src.rows),
                   openvx_equalize_hist(src, dst))
    
        Mutex histogramLockInstance;
    
        const int hist_sz = EqualizeHistCalcHist_Invoker::HIST_SZ;
        int hist[hist_sz] = {0,}; // 直方图表
        int lut[hist_sz]; // 查找表
    
        EqualizeHistCalcHist_Invoker calcBody(src, hist, &histogramLockInstance); // 建立直方图操作
        EqualizeHistLut_Invoker      lutBody(src, dst, lut); // 查找表操作
        cv::Range heightRange(0, src.rows); // 总行数的范围
    
        if(EqualizeHistCalcHist_Invoker::isWorthParallel(src))
            parallel_for_(heightRange, calcBody);
        else
            calcBody(heightRange);
    
        int i = 0;
        while (!hist[i]) ++i; // 跳过那些像素个数为0的
    
        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); // 前面的像素在直方图横坐标上能占多宽，注意这里是uchar的整数
        }
    	
    	// 执行均衡化操作
        if(EqualizeHistLut_Invoker::isWorthParallel(src))
            parallel_for_(heightRange, lutBody);
        else
            lutBody(heightRange);
    }

可以看出，equalizeHist 定义了EqualizeHistCalcHist\_Invoker calcBody(src, hist, &histogramLockInstance)来建立直方图，同时还定义了EqualizeHistLut\_Invoker  lutBody(src, dst, lut)来从查找表构建目标图像。这两个类都继承了抽象类cv::ParallelLoopBody，该类是一个专门为并行处理设计的类，其对象会在parallel\_for\_impl中传递给Concurrency::parallel\_for函数。并行处理能更充分地发挥设备的硬件性能，所以速度非常快。

因为涉及到并行处理，我们先理解一下OpenCV并行处理的过程，同时也可以借鉴一下Intel公司在并行处理上的编码风格和艺术。

并行处理的真正实现函数是parallel\_for\_impl，其中定义的ParallelLoopBodyWrapperContext在始化时起到联接上下文的作用，在并行运行时本身没有实质内容。

在看源码之前，我们先理一下运行的过程和顺序，equilizeHist有两处用到了parallel\_for\_

EqualizeHistCalcHist_Invoker calcBody(src, hist, &histogramLockInstance);
    EqualizeHistLut_Invoker      lutBody(src, dst, lut);
    …
    if(EqualizeHistCalcHist_Invoker::isWorthParallel(src))
            parallel_for_(heightRange, calcBody);
        else
            calcBody(heightRange);
    …
        if(EqualizeHistLut_Invoker::isWorthParallel(src))
            parallel_for_(heightRange, lutBody);
        else
            lutBody(heightRange);
    …

前面讲过，calcBody负责计算直方图，lutBody负责填充目标图像，这两个parallel\_for\_执行的过程是完全一样的。在parallel\_for\_impl函数中，并行运行的系统API是

Concurrency::parallel\_for(stripeRange.start, stripeRange.end, pbody);

设置的处理函数是ProxyLoopBody pbody，这两个Body执行的路线（关系）就是ProxyLoopBody pbody(ctx)----> ParallelLoopBodyWrapperContext ctx(body)----> Body, 其实现都在相关类的()重载函数中。

首先pbody运行时，实例ProxyLoopBody pbody(ctx)被执行，即

class ProxyLoopBody : public ParallelLoopBodyWrapper
        {
            void operator ()(int i) const
            {
                this->ParallelLoopBodyWrapper::operator()(cv::Range(i, i + 1));
            }
        };

其中range就是

Concurrency::parallel\_for(stripeRange.start, stripeRange.end, pbody);

传入的参数范围，系统会根据这个范围逐个处理，如果启用TBB，就是每次处理一个小范围，直到全部处理完毕。

而后，这个operator()重载函数执行时，其被继承的ParallelLoopBodyWrapper类的()重载函数会被调用，

class ParallelLoopBodyWrapper : public cv::ParallelLoopBody
        {  
            void operator()(const cv::Range& sr) const
           {
                // propagate main thread state
                cv::theRNG() = ctx.rng; // 这个就是ParallelLoopBodyWrapperContext里初始化时得到的线程状态
    
                cv::Range r;
                cv::Range wholeRange = ctx.wholeRange;
                int nstripes = ctx.nstripes;
                r.start = (int)(wholeRange.start +
                                ((uint64)sr.start*(wholeRange.end - wholeRange.start) + nstripes/2)/nstripes);
                r.end = sr.end >= nstripes ? wholeRange.end : (int)(wholeRange.start +
                                ((uint64)sr.end*(wholeRange.end - wholeRange.start) + nstripes/2)/nstripes);
    
                (*ctx.body)(r); // 计算直方图 或 填充目标图像
    
                if (!ctx.is_rng_used && !(cv::theRNG() == ctx.rng))
                    ctx.is_rng_used = true;
            }
    }

其中，在构建直方图时，ctx.body就是EqualizeHistCalcHist\_Invoker的实例，

(\*ctx.body)(r);

会执行下面这段直方图的计算

class EqualizeHistCalcHist_Invoker : public cv::ParallelLoopBody
    {
        void operator()( const cv::Range& rowRange ) const
        {。。。}
    }

在构建目标图像时，ctx.body就是EqualizeHistLut\_Invoker的实例，(\*ctx.body)(r)就相当于运行，

class EqualizeHistLut_Invoker : public cv::ParallelLoopBody
    {
        void operator()( const cv::Range& rowRange ) const
        {
            。。。
        }
    }

下面是parallel\_for\_impl的源码（删除了一些与本文主旨无关的），很简单，

static void parallel_for_impl(const cv::Range& range, const cv::ParallelLoopBody& body, double nstripes)
    {
        if ((numThreads < 0 || numThreads > 1) && range.end - range.start > 1)
        {
            ParallelLoopBodyWrapperContext ctx(body, range, nstripes); // 这里会初始化ctx->nstripes
            ProxyLoopBody pbody(ctx); 
            cv::Range stripeRange = pbody.stripeRange();     // 就是上面ctx初始化的那个ctx->nstripes
            if( stripeRange.end - stripeRange.start == 1 )
            {
                body(range);
                return;
            }
    
            if(!pplScheduler || pplScheduler->Id() == Concurrency::CurrentScheduler::Id())
            {
                Concurrency::parallel_for(stripeRange.start, stripeRange.end, pbody); // 打开多线程
            }
            else
            {
                pplScheduler->Attach();
                Concurrency::parallel_for(stripeRange.start, stripeRange.end, pbody);
                Concurrency::CurrentScheduler::Detach();
            }
        }
    }

## OpenCV均衡化的具体算法详解 ##

在具体算法上，用一个实在的例子来讲比较容易理解。

OpenCV是这样处理的，

先建立直方图hist\[256\]，比如某图片的灰度图，

Hist\[0~15\] = 0，表示这几个颜色没有用到

hist\[16\] = 1, 表示颜色为16的有1个像素，

hist\[17\] = 2, 表示颜色为17的有10个像素，

hist\[18\] = 5, 表示颜色为18的有29个像素，

。。。

然后，用equalizeHist计算scale,这个scale的物理意义，大致可以理解为每个像素在直方图上的横坐标上能占多宽。

比如得到scale=0.1后建立查找表lut\[256\]，就是这样的，

Lut\[0~15\] = 0

Lut\[16\] = (int)0.1\*1 = 0

Lut\[17\] = (int)0.1\*(1+10)=1

Lut\[18\] = (int)0.1\*(1+10+29) = 4

。。。

再然后，在新的图像中，

原颜色为16的 改成 新色 0

原颜色为17的 改成 新色1

原颜色为18的 改成 新色4

。。。

这样，灰度值的动态范围得到了扩大，对比度得到加强，图像可以进一步处理。

## 源码 ##

### EqualizeHistCalcHist\_Invoker ###

作用：进行直方图统计

class EqualizeHistCalcHist_Invoker : public cv::ParallelLoopBody
    {
    public:
        enum {HIST_SZ = 256};
    
        EqualizeHistCalcHist_Invoker(cv::Mat& src, int* histogram, cv::Mutex* histogramLock)
            : src_(src), globalHistogram_(histogram), histogramLock_(histogramLock)
        { }
    
        void operator()( const cv::Range& rowRange ) const 
        {
            int localHistogram[HIST_SZ] = {0, }; 
    
            const size_t sstep = src_.step;
    
            int width = src_.cols;
            int height = rowRange.end - rowRange.start;
    
            if (src_.isContinuous())
            {
                width *= height;
                height = 1;
            }
    
            for (const uchar* ptr = src_.ptr<uchar>(rowRange.start); height--; ptr += sstep)
            {
                int x = 0;
    		   // 下面这个for循环，就是统计图片的直方图
    		   // 每次统计4个点，
                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)
                    localHistogram[ptr[x]]++;
            }
    
            cv::AutoLock lock(*histogramLock_);
    
    		// 把结果放入globalHistogram_
            for( int i = 0; i < HIST_SZ; i++ )
                globalHistogram_[i] += localHistogram[i];
        }
    
        static bool isWorthParallel( const cv::Mat& src )
        {
            return ( src.total() >= 640*480 );
        }
    
    private:
        EqualizeHistCalcHist_Invoker& operator=(const EqualizeHistCalcHist_Invoker&);
    
        cv::Mat& src_;
        int* globalHistogram_;
        cv::Mutex* histogramLock_;
    };

### EqualizeHistLut\_Invoker ###

作用：执行直方图的均衡化操作

class EqualizeHistLut_Invoker : public cv::ParallelLoopBody
    {
    public:
        EqualizeHistLut_Invoker( cv::Mat& src, cv::Mat& dst, int* lut )
            : src_(src),
              dst_(dst),
              lut_(lut)
        { }
    
        void operator()( const cv::Range& rowRange ) const
        {
            const size_t sstep = src_.step;
            const size_t dstep = dst_.step;
    
            int width = src_.cols;
            int height = rowRange.end - rowRange.start;
            int* lut = lut_;
    
            if (src_.isContinuous() && dst_.isContinuous())
            {
                width *= height;
                height = 1;
            }
    
            const uchar* sptr = src_.ptr<uchar>(rowRange.start);
            uchar* dptr = dst_.ptr<uchar>(rowRange.start);
    	
            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 )
        {
            return ( src.total() >= 640*480 );
        }
    
    private:
        EqualizeHistLut_Invoker& operator=(const EqualizeHistLut_Invoker&);
    
        cv::Mat& src_;
        cv::Mat& dst_;
        int* lut_;
    };

解析完毕

[70]: https://img-blog.csdn.net/20180815115917356?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3Rhbm14MjE5/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70
[70 1]: /images/20220516/3f873b918a6d4a349d10e98d970a1945.png
[70 2]: https://img-blog.csdn.net/20180815120010121?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3Rhbm14MjE5/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70
[70 3]: /images/20220516/7fededc91b044b0ab95bfe2975fe675b.png
[70 4]: https://img-blog.csdn.net/20180815120027490?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3Rhbm14MjE5/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70
[70 5]: /images/20220516/2ded8100dde34e7a86fbf800b5888822.png