直方图均衡化
之前的文章中陆续介绍了OpenCV的编译，色彩空间以及滤波器，甚至DNN的简单介绍，挖了不少坑，目前很多都还没有填上，东西很多，也很杂乱。为了方便读者学习，从本文开始，我将从OpenCV的基本的图像处理算法开始，逐步系统的介绍OpenCV的各个模块的功能。本文先从直方图均衡化开始介绍。网上关于OpenCV API使用方法的文章非常多，但是对于背后的算法原理介绍就比较少了，所以在后续的文章中将侧重于OpenCV算法原理以及OpenCV算法实现源码的介绍。
直方图均衡化是图像处理中非常基本的算法，但是却有着非常重要的作用，它是多种空间域处理技术的基础，除了可以提供图像的统计信息外，还可用于图像增强，图像压缩，图像分割等等，由于其算法软件计算简单，且有助于商业硬件实现，因此已经成为实时图像处理的流行工具。
算法原理简介
直方图均衡化是一种通过拉伸像素强度分布来增强图像对比度的方法，直方图均衡化的算法原理如名字一样，先求直方图，而后做均衡化，分为两部分。所以算法原理方面分为两小节。
直方图
关于直方图函数的推来过程在数字图像处理的教材上有着详细的推到，这里不做赘述。这里直接给出计算公式：
其中，MN是图像的像素总数量，nk表示灰度为rk的像素个数，L表示灰度级的数量，例如8bit图像，灰度级为256；与rk相对的pr(rk)图形通常称为直方图。
也就是说，图像的直方图就是图像中每个灰度级的像素在图像中出现的频率。获得直方图以后，就需要进行均衡化处理。
均衡化
图像均衡化处理其实是积分的过程，对于数字图像的离散值便是求和。计算公式如下：
以上是均衡化的计算公式。其实就是(L-1)与像素出现频率累加和的乘积。
这个算法很简单，以上就是算法的基本计算原理。接下我们看看OpenCV的源码实现方式。
OpenCV源码
虽然算法很简单，但是在OpenCV工程源代码代码中也并不简单，其中涉及到各个平台的优化等等，在之前的文章中我们也介绍过一些，和以前一样，我们不深入分析各个平台优化的代码，只针对最通用的c++实现给出分析，目的也在于能够更好的理解算法的实现。关于优化会在后续高性能计算部分介绍。
API介绍
void cv::equalizeHist(InputArray src,
OutputArray dst)相当的简单，直接给出输入输出即可。接下来看API源码。
源码
path： opencv/modules/imgproc/src/histogram.cpp +3330
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;
// HIST_SZ = 256;
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;
//查找图片中第一个频率不为0的灰度级i
while (!hist[i]) ++i;
int total = (int)src.total();
//如果第一个出现频率不为0的灰度级的出现的频次与输入图像的像素总数相等；
//也就是说整幅图片就一个灰度级，那么均衡化的结果为dst所有像素为i
if (hist[i] == total)
{
dst.setTo(i);
return;
}
//这里就是上文均衡化公式的(L-1)/MN
//需要注意这里分母减去hist[i];这里是为了剔除前i个连续的频次为0的灰度级
//大家可以自行推到，前n个连续的频次为0的点在累加过程中是没有作用的，所以剔除
float scale = (hist_sz - 1.f)/(total - hist[i]);
int sum = 0;
//计算Lut
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);
}整个程序的计算流程我在主程序中做了详细的注释。这里对于关键点进行重点介绍：
我们可以看到程序的主要计算分为3个部分：
- 计算直方图(calcBody)
- 计算查找表：查找表是一个数组，索引就是灰度级，即0-255，对应的值就是该灰度级均衡化以后的灰度级。例如Lut[24]=30；表示图片中像素值为24的像素，均衡化以后，像素值为30；因此计算查找的过程其实就是我们上文公式中均衡化的过程。
- 均衡化，这里就是运用查找表，讲原图像中的像素替换为均衡化以后的像素，写如dst；
以上便是OpenCV中直方图均衡化实现的细节。下面展开介绍直方图的计算和均衡化两个过程，查找表的计算使用的就是算法原理中均衡化的公式，在主程序中进行了详细注释，不再赘述。
直方图计算源码：
path： opencv/modules/imgproc/src/histogram.cpp +3114
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 CV_OVERRIDE
{
int localHistogram[HIST_SZ] = {0, };
const size_t sstep = src_.step;
int width = src_.cols;
int height = rowRange.end - rowRange.start;
//如果内存是连续的，就把二维空间转换为1维
if (src_.isContinuous())
{
width *= height;
height = 1;
}
//for循环遍历src，计数各个灰度级出现的频次
for (const uchar* ptr = src_.ptr<uchar>(rowRange.start); height--; ptr += sstep)
{
int x = 0;
//在width方西上对for循环进行4阶展开
//如果内存连续，则height方向只循环一次
for (; x <= width - 4; x += 4)
{
//读取原图像像素值
int t0 = ptr[x], t1 = ptr[x+1];
//对应灰度级频次计数加以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_);
//汇总每个线程计算分块的局部直方图，变为整个图像的全局直方图
for( int i = 0; i < HIST_SZ; i++ )
globalHistogram_[i] += localHistogram[i];
}
//该函数用于判断是否使用多线程或者说并行for循环
static bool isWorthParallel( const cv::Mat& src )
{
//如果图像大于640*480就使用并行for循环
return ( src.total() >= 640*480 );
}
private:
EqualizeHistCalcHist_Invoker& operator=(const EqualizeHistCalcHist_Invoker&);
cv::Mat& src_;
int* globalHistogram_;
cv::Mutex* histogramLock_;
};以上是对于图像直方图计算函数的详细注释。需要留意的是，OpenCV除了使用并行for循环这样的优化以外，还是用了循环展开这样的指令级向量优化。
均衡化源码：
path： opencv/modules/imgproc/src/histogram.cpp +3114
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 CV_OVERRIDE
{
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_;
//内存连续转换为1维空间
if (src_.isContinuous() && dst_.isContinuous())
{
width *= height;
height = 1;
}
const uchar* sptr = src_.ptr<uchar>(rowRange.start);
uchar* dptr = dst_.ptr<uchar>(rowRange.start);
//height方向外循环，若内存连续，height方向循环一次
for (; height--; sptr += sstep, dptr += dstep)
{
int x = 0;
//width方向内循环，循环4阶展开
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]];
}
}
//判断是否使用并行for循环
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_;
};以上是均衡化过程的源代码，进行了详细的注释。大家可以自行阅读。
整个OpenCV的直方图均衡化源码分析就到此结束了。下一篇文章分析OpenCV的空域滤波，可能需要好几篇分析，毕竟空域滤波的内容很多。关于API使用，由于不是本系列文章的重点，所以会在空域滤波结束后，专门写一片文章，统一展示多个API的使用。
经过本片文章分析，可以至少学到使用循环展开来加速自己的代码。不同机器体处理器架构不同，可能展开的阶次略有区别。过渡展开和展开不足都会影响程序性能。
欢迎关注公众号：计算机视觉与高性能计算(to_konw)
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
