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

理论基础

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

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

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

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()  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(rowRange.start); height--; ptr += sstep)
        {
            int x = 0;
		   // 下面这个for循环，就是统计图片的直方图
		   // 每次统计4个点，
            for (; x = 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(rowRange.start);
        uchar* dptr = dst_.ptr(rowRange.start);
	
        for (; height--; sptr += sstep, dptr += dstep)
        {
            int x = 0;
            for (; x = 640*480 );
    }

private:
    EqualizeHistLut_Invoker& operator=(const EqualizeHistLut_Invoker&);

    cv::Mat& src_;
    cv::Mat& dst_;
    int* lut_;
};
 
解析完毕