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() <= 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_;
};
解析完毕
