图像直方图及均衡方法总结（二）自适应直方图均衡AHE,CLAHE

在http://blog.csdn.net/piaoxuezhong/article/details/78269439一文中对图像直方图及直方图均衡做了总结，由于篇幅原因，后面的自适应直方图均衡部分单独在本篇总结一下。

自适应直方图均衡化(Adaptive histgram equalization/AHE)

按照维基百科里的说法，“It differs from ordinary histogram equalization in the respect that the adaptive method computes several
histograms, each corresponding to a distinct section of the image, and uses them to redistribute the lightness values of the image. It is therefore suitable for improving the local contrast and enhancing the definitions of edges in each region of an image.”

也就是说：AHE算法与经典算法的不同点在于它通过计算图像多个局部区域的直方图，并重新分布亮度，以此改变图像对比度。所以，该算法更适合于提高图像的局部对比度和细节部分。不过呢，AHE存在过度放大图像中相对均匀区域的噪音的问题。

直方图均衡经典算法对整幅图像的像素使用相同的变换，对于像素值分布比较均衡的图像来说，经典算法的效果不错。但是如果图像中包括明显亮的或者暗的区域，在这些部分的对比度并不能得到增强。

AHE算法通过对局部区域的直方图均衡解决上述问题。最简单的形式就是每个像素通过对其周边一个矩形范围内像素的直方图进行均衡化，均衡的方式使用普通的均衡化算法，变换函数同像素周边的累积分布函数（CDF）成比例。当然局部直方图还可以采取其他方式，局部直方图处理大致有3种实现方法：

1）将原始图片划分成不重叠的子块，在每个子块内做直方图处理，该方法简单但输出图像会有块效应；

2）类似模板卷积的方式，以待处理的点为中心，取其邻域为子块，在子块内做直方图处理，处理结果仅映射到该点可以消除块效应，但需要对每个点计算一次直方图处理，效率低；

3）前两种方法的结合版，不再逐像素移动，步长小于子块宽度以确保两个相邻子块有重叠；每个子块做直方图映射后的结果赋值给子块内所有点，这样每个点会有多次赋值，最终的取值为这些赋值的均值；

由于AHE算法存在放大噪声等问题，后来的研究者们发明了限制对比度的自适应直方图均衡算法（CLAHE），关于CLAHE的说明和使用可以找到很多，我这里尽量全面的汇总一下。

限制对比度自适应直方图均衡（Contrast Limited Adaptive histgram equalization/CLAHE)

CLAHE在局部直方图均衡化（又称自适应直方图均衡化AHE）的基础上，对每个子块直方图做了限制，很好的控制了AHE带来的噪声。（后面对维基百科里的部分解释搬砖...）

CLAHE与AHE的不同主要在于其对于对比度的限幅，在CLAHE中，对于每个小区域都必须使用对比度限幅，用来克服AHE的过度放大噪音的问题。 

算法通过限制AHE的对比增强程度来达到效果。对于一个像素值，其周边的对比度放大主要是由变换函数的斜吕决定的，该斜度与区域累积分布函数的斜度成比例。在计算CDF前，CLAHE通过用预先定义的阈值来裁剪直方图以达到限制放大倍数的目的。该算法的优势在于它不是选择直接忽略掉那些超出限幅的部分，而是将这些裁剪掉的部分均匀分布到直方图的其他部分。



为了提升计算速度及去除分块处理导致的块边缘过渡不平衡等问题，可以使用插值方法。关于算法的介绍和描述，请参见文章。对于算法实现，可以选择：

（1）维基百科里给出了CLAHE的C语言版本的实现：http://www.realtimerendering.com/resources/GraphicsGems/gemsiv/clahe.c
（2）matlab给出的函数adapthisteq，想看源码的可以在matlab里ctrl+D，这里不再附了。这里举个实例：

clear all,close all,clc;
img = imread('../image/7.jpg');
figure,
subplot(1,2,1),imshow(img);title('原图');
if length(size(img))>2
rimg = img(:,:,1);
gimg = img(:,:,2);
bimg = img(:,:,3);
resultr = adapthisteq(rimg);
resultg = adapthisteq(gimg);
resultb = adapthisteq(bimg);
result = cat(3, resultr, resultg, resultb);
subplot(1,2,2),imshow(result);title('CLAHE处理后');
end

当然，还有建议在Lab空间中操作的，具体实现就不再写了，空间转换一下就好了。

（3）opencv里也有CLAHE的实现，目录为：..\opencv\sources\modules\imgproc\src\clahe.cpp

我这里附下源码：http://blog.csdn.net/panda1234lee/article/details/52852765

// CLAHE

namespace
{
class CLAHE_CalcLut_Body : public cv::ParallelLoopBody
{
public:
CLAHE_CalcLut_Body(const cv::Mat& src, cv::Mat& lut, cv::Size tileSize, int tilesX, int clipLimit, float lutScale) :
src_(src), lut_(lut), tileSize_(tileSize), tilesX_(tilesX), clipLimit_(clipLimit), lutScale_(lutScale)
{
}

void operator ()(const cv::Range& range) const;

private:
cv::Mat src_;
mutable cv::Mat lut_;

cv::Size tileSize_;
int tilesX_;
int clipLimit_;
float lutScale_;
};

// 计算直方图查找表
void CLAHE_CalcLut_Body::operator ()(const cv::Range& range) const
{
const int histSize = 256;

uchar* tileLut = lut_.ptr(range.start);
const size_t lut_step = lut_.step;  // size = tilesX_*tilesY_ * lut_step

// Range(0, tilesX_ * tilesY_)，全图被分为tilesX_*tiles_Y个块
for (int k = range.start; k < range.end; ++k, tileLut += lut_step)
{
// (tx, ty)表示当前所在是哪一块
// (0, 0) (1, 0)...(tilesX_-1, 0)
// (0, 1) (1, 1)...(tilesX_-1, 1)
// ...
// (0, tilesY_-1)... (tilesX_-1, tilesY_-1)
const int ty = k / tilesX_;
const int tx = k % tilesX_;

// retrieve tile submatrix
// 注意：tileSize.width表示分块的宽度，tileSize.height表示分块高度
cv::Rect tileROI;
tileROI.x = tx * tileSize_.width;   // 换算为全局坐标
tileROI.y = ty * tileSize_.height;
tileROI.width = tileSize_.width;
tileROI.height = tileSize_.height;

const cv::Mat tile = src_(tileROI);

// calc histogram
int tileHist[histSize] = { 0, };
// 统计 ROI 的直方图
int height = tileROI.height;
const size_t sstep = tile.step;
for (const uchar* ptr = tile.ptr<uchar>(0); height--; ptr += sstep)
{
int x = 0;
for (; x <= tileROI.width - 4; x += 4)
{
int t0 = ptr[x], t1 = ptr[x + 1];
tileHist[t0]++; tileHist[t1]++;
t0 = ptr[x + 2]; t1 = ptr[x + 3];
tileHist[t0]++; tileHist[t1]++;
}

for (; x < tileROI.width; ++x)
tileHist[ptr[x]]++;
}

// clip histogram
if (clipLimit_ > 0)
{
// how many pixels were clipped
int clipped = 0;
for (int i = 0; i < histSize; ++i)
{
// 超过裁剪阈值
if (tileHist[i] > clipLimit_)
{
clipped += tileHist[i] - clipLimit_;
tileHist[i] = clipLimit_;
}
}

// redistribute clipped pixels
int redistBatch = clipped / histSize;
int residual = clipped - redistBatch * histSize;

// 平均分配裁剪的差值到所有直方图
for (int i = 0; i < histSize; ++i)
tileHist[i] += redistBatch;

// 处理差值的余数
for (int i = 0; i < residual; ++i)
tileHist[i]++;
}

// calc Lut
int sum = 0;
for (int i = 0; i < histSize; ++i)
{
// 累加直方图
sum += tileHist[i];
tileLut[i] = cv::saturate_cast<uchar>(sum * lutScale_);   // static_cast<float>(histSize - 1) / tileSizeTotal
}
}
}

class CLAHE_Interpolation_Body : public cv::ParallelLoopBody
{
public:
CLAHE_Interpolation_Body(const cv::Mat& src, cv::Mat& dst, const cv::Mat& lut, cv::Size tileSize, int tilesX, int tilesY) :
src_(src), dst_(dst), lut_(lut), tileSize_(tileSize), tilesX_(tilesX), tilesY_(tilesY)
{
}

void operator ()(const cv::Range& range) const;

private:
cv::Mat src_;
mutable cv::Mat dst_;
cv::Mat lut_;

cv::Size tileSize_;
int tilesX_;
int tilesY_;
};

// 根据相邻4块的直方图插值
void CLAHE_Interpolation_Body::operator ()(const cv::Range& range) const
{
const size_t lut_step = lut_.step;
// Range(0, src.rows)
for (int y = range.start; y < range.end; ++y)
{
const uchar* srcRow = src_.ptr<uchar>(y);
uchar* dstRow = dst_.ptr<uchar>(y);

const float tyf = (static_cast<float>(y) / tileSize_.height) - 0.5f;

int ty1 = cvFloor(tyf);
int ty2 = ty1 + 1;

// 差值作为插值的比例
const float ya = tyf - ty1;

ty1 = std::max(ty1, 0);
ty2 = std::min(ty2, tilesY_ - 1);

const uchar* lutPlane1 = lut_.ptr(ty1 * tilesX_);   // 当前块的直方图
const uchar* lutPlane2 = lut_.ptr(ty2 * tilesX_);   // 向下一块的直方图

for (int x = 0; x < src_.cols; ++x)
{
const float txf = (static_cast<float>(x) / tileSize_.width) - 0.5f;

int tx1 = cvFloor(txf);
int tx2 = tx1 + 1;

// 差值作为插值的比例
const float xa = txf - tx1;

tx1 = std::max(tx1, 0);
tx2 = std::min(tx2, tilesX_ - 1);

// src_.ptr<uchar>(y)[x]
const int srcVal = srcRow[x];

// 索引 LUT
const size_t ind1 = tx1 * lut_step + srcVal;
const size_t ind2 = tx2 * lut_step + srcVal;  // 向右一块的直方图

float res = 0;
// 根据直方图的值进行插值
// lut_.ptr(ty1 * tilesX_)[tx1 * lut_step + srcVa] => lut_[ty1][tx1][srcVal]
res += lutPlane1[ind1] * ((1.0f - xa) * (1.0f - ya));
res += lutPlane1[ind2] * ((xa) * (1.0f - ya));
res += lutPlane2[ind1] * ((1.0f - xa) * (ya));
res += lutPlane2[ind2] * ((xa) * (ya));

dstRow[x] = cv::saturate_cast<uchar>(res);
}
}
}

class CLAHE_Impl : public cv::CLAHE
{
public:
CLAHE_Impl(double clipLimit = 40.0, int tilesX = 8, int tilesY = 8);

cv::AlgorithmInfo* info() const; // Algorithm类工厂方法封装相关

void apply(cv::InputArray src, cv::OutputArray dst);

void setClipLimit(double clipLimit);
double getClipLimit() const;

void setTilesGridSize(cv::Size tileGridSize);
cv::Size getTilesGridSize() const;

void collectGarbage();

private:
double clipLimit_;
int tilesX_;
int tilesY_;

cv::Mat srcExt_;
cv::Mat lut_;
};

CLAHE_Impl::CLAHE_Impl(double clipLimit, int tilesX, int tilesY) :
clipLimit_(clipLimit), tilesX_(tilesX), tilesY_(tilesY)
{
}
// Algorithm类工厂方法封装相关
//CV_INIT_ALGORITHM(CLAHE_Impl, "CLAHE",
//  obj.info()->addParam(obj, "clipLimit", obj.clipLimit_);
//obj.info()->addParam(obj, "tilesX", obj.tilesX_);
//obj.info()->addParam(obj, "tilesY", obj.tilesY_))

void CLAHE_Impl::apply(cv::InputArray _src, cv::OutputArray _dst)
{
cv::Mat src = _src.getMat();

CV_Assert(src.type() == CV_8UC1);

_dst.create(src.size(), src.type());
cv::Mat dst = _dst.getMat();

const int histSize = 256;
// 准备 LUT，tilesX_*tilesY_个块，每个块都有256个柱子的直方图
lut_.create(tilesX_ * tilesY_, histSize, CV_8UC1);

cv::Size tileSize;
cv::Mat srcForLut;

// 如果分块刚好（整除）
if (src.cols % tilesX_ == 0 && src.rows % tilesY_ == 0)
{
tileSize = cv::Size(src.cols / tilesX_, src.rows / tilesY_);
srcForLut = src;
}
// 否则对原图进行扩充
else
{
cv::copyMakeBorder(src, srcExt_, 0, tilesY_ - (src.rows % tilesY_), 0, tilesX_ - (src.cols % tilesX_), cv::BORDER_REFLECT_101);

tileSize = cv::Size(srcExt_.cols / tilesX_, srcExt_.rows / tilesY_);
srcForLut = srcExt_;
}

const int tileSizeTotal = tileSize.area();
const float lutScale = static_cast<float>(histSize - 1) / tileSizeTotal;      // △

// 计算实际的clipLimit
int clipLimit = 0;
if (clipLimit_ > 0.0)
{
clipLimit = static_cast<int>(clipLimit_ * tileSizeTotal / histSize);
clipLimit = std::max(clipLimit, 1);
}

// 分块并行计算： LUT
CLAHE_CalcLut_Body calcLutBody(srcForLut, lut_, tileSize, tilesX_, clipLimit, lutScale);
cv::parallel_for_(cv::Range(0, tilesX_ * tilesY_), calcLutBody);

// 分块并行计算： 根据直方图插值
CLAHE_Interpolation_Body interpolationBody(src, dst, lut_, tileSize, tilesX_, tilesY_);
cv::parallel_for_(cv::Range(0, src.rows), interpolationBody);
}

void CLAHE_Impl::setClipLimit(double clipLimit)
{
clipLimit_ = clipLimit;
}

double CLAHE_Impl::getClipLimit() const
{
return clipLimit_;
}

void CLAHE_Impl::setTilesGridSize(cv::Size tileGridSize)
{
tilesX_ = tileGridSize.width;
tilesY_ = tileGridSize.height;
}

cv::Size CLAHE_Impl::getTilesGridSize() const
{
return cv::Size(tilesX_, tilesY_);
}

void CLAHE_Impl::collectGarbage()
{
srcExt_.release();
lut_.release();
}
}

cv::Ptr<cv::CLAHE> cv::createCLAHE(double clipLimit, cv::Size tileGridSize)
{
return new CLAHE_Impl(clipLimit, tileGridSize.width, tileGridSize.height);
}

opencv实现CLAHE的主要步骤，可以分为如下几步:

step1. 扩展图像边界，使其能够刚好切分为若干子块，假设每个子块面积为tileSizeTotal，子块系数 lutScale = 255.0 / tileSizeTotal，对预设的limit做处理：limit = MAX(1, limit * tileSizeTotal / 256)；

step2. 对每个子块，计算直方图；

step3. 对每个子块直方图的每个灰度级，使用预设的limit值做限定，同时统计整个直方图超过limit的像素数；

step4. 计算每个子块的lut累积直方图tileLut，tileLut[i] = sum[i] * lutScale，sum[i] 是累积直方图，lutScale确保tileLut取值在[0, 255]；

step5. 遍历原始图像每个点，考虑该点所在子块及右、下、右下一共4个子块的tileLut，以原始灰度值为索引得到4个值，然后做双线性插值得该点变换后的灰度值；
下面是我使用opencv对彩色图像的CLAHE实现，注意程序里将图像从RGB空间转为Lab空间：
#include "opencv2/opencv.hpp"
#include <iostream>

using namespace cv;
using namespace std;

void main()
{
Mat inp_img = cv::imread("D:\\proMatlab\\vessel_edge_extration\\image\\7.jpg");
if (!inp_img.data) {
cout << "Something Wrong";
}
namedWindow("Input Image", CV_WINDOW_AUTOSIZE);
imshow("Input Image", inp_img);

Mat clahe_img;
cvtColor(inp_img, clahe_img, CV_BGR2Lab);
vector<cv::Mat> channels(3);
split(clahe_img, channels);

Ptr<cv::CLAHE> clahe = createCLAHE();
// 直方图的柱子高度大于计算后的ClipLimit的部分被裁剪掉，然后将其平均分配给整张直方图
clahe->setClipLimit(4.); // (int)(4.*(8*8)/256)
//clahe->setTilesGridSize(Size(8, 8)); // 将图像分为8*8块
Mat dst;
clahe->apply(channels[0], dst);
dst.copyTo(channels[0]);
//
clahe->apply(channels[1], dst);
dst.copyTo(channels[1]);
clahe->apply(channels[2], dst);
dst.copyTo(channels[2]);
merge(channels, clahe_img);

Mat image_clahe;
cvtColor(clahe_img, image_clahe, CV_Lab2BGR);

namedWindow("CLAHE Image", CV_WINDOW_AUTOSIZE);
imshow("CLAHE Image", image_clahe);
imwrite("out.jpg", image_clahe);
waitKey();
system("pause");
}原始图同在matlab里处理的一样，这里只附处理后的结果：

CLAHE算法的应用

目前，CLAHE算法在图像处理领域有较多的应用，典型的有图像去雾，水下图像处理，低照度图像处理等方面。当然具体到实际应用，算法是需要调参和优化的。

参考

https://en.wikipedia.org/wiki/Adaptive_histogram_equalization#Contrast_Limited_AHE
http://blog.csdn.net/u010839382/article/details/49584181 http://blog.csdn.net/fightingforcv/article/details/47661463 http://blog.csdn.net/linear_luo/article/details/52527563 http://blog.csdn.net/baimafujinji/article/details/50660189 http://www.cnblogs.com/Imageshop/archive/2013/04/07/3006334.html