直方图是图像色彩统计特征的抽象表述。基于直方图可以实现很多有趣的算法。例如，图像增强中利用直方图来调整图像的对比度、有人利用直方图来进行大规模无损数据隐藏、还有人利用梯度直方图HOG来构建图像特征进而实现目标检测。本节我们就来讨论重要的直方图均衡化算法，说它重要是因为以此为基础后续又衍生出了许多实用而有趣的算法。
Histogram equalization
如果一幅图像的像素灰度值在一个过于有限的范围内聚集，那么图像的程序效果即会很糟糕，直接观感就是对比度很弱。下图来自维基百科，第一幅图的直方图分布非常不均衡。如果把直方图均匀地延展到整个分布域内，则图像的效果显得好了很多。
Matlab中提供了现成的函数“histeq()”来实现图像的直方图均衡。但为了演示说明算法的原理，下面我将在Matlab中自行编码实现图像的直方图均衡。通过代码来演示这个算法显然更加直观，更加易懂。当然，其实我还不得不感叹，如果仅仅是作为图像算法研究之用，Matlab确实非常好用。
首先读入图像，并将其转化为灰度图。然后提取图像的长和宽。
image = imread('Unequalized_Hawkes_Bay_NZ.jpg');
Img = rgb2gray(image);
[height,width]=size(image);然后绘制一下原始图像的直方图。
[counts1, x] = imhist(Img,256);
counts2 = counts1/height/width;
stem(x, counts2);
统计每个灰度的像素值累计数目。
NumPixel = zeros(1,256);%统计各灰度数目，共256个灰度级
for i = 1:height
for j = 1: width
%对应灰度值像素点数量增加一
%因为NumPixel的下标是从1开始，但是图像像素的取值范围是0~255，所以用NumPixel(Img(i,j) + 1)
NumPixel(Img(i,j) + 1) = NumPixel(Img(i,j) + 1) + 1;
end
end然后将频数值算为频率
ProbPixel = zeros(1,256);
for i = 1:256
ProbPixel(i) = NumPixel(i) / (height * width * 1.0);
end再用函数cumsum来计算cdf，并将频率（取值范围是0.0~1.0）映射到0~255的无符号整数。
CumuPixel = cumsum(ProbPixel);
CumuPixel = uint8(255 .* CumuPixel + 0.5);
直方图均衡。赋值语句右端，Img(i,j)被用来作为CumuPixel的索引。比如Img(i,j) = 120，则从CumuPixel中取出第120个值作为Img(i,j) 的新像素值。
for i = 1:height
for j = 1: width
Img(i,j) = CumuPixel(Img(i,j));
end
end
最后显示新图像的直方图。
imshow(Img);
[counts1, x] = imhist(Img,256);
counts2 = counts1/height/width;
stem(x, counts2);
当然，上述讨论的是灰度图像的直方图均衡。对于彩色图像而言，你可能会想到分别对R、G、B三个分量来做处理，这也确实是一种方法。但有些时候，这样做很有可能导致结果图像色彩失真。因此有人建议将RGB空间转换为HSV之后，对V分量进行直方图均衡处理以保存图像色彩不失真。下面我们来做一些对比实验。待处理图像是标准的图像处理测试用图couple图，如下所示。
首先，我们分别处理R、G、B三个分量，为了简便我们直接使用matlab中的函数histeq()。
a = imread('couple.tiff');
R = a(:,:,1);
G = a(:,:,2);
B = a(:,:,3);
R = histeq(R, 256);
G = histeq(G, 256);
B = histeq(B, 256);
a(:,:,1) = R;
a(:,:,2) = G;
a(:,:,3) = B;
imshow(a)下面的代码使用了另外一种方式，即将色彩空间转换到HSV后，对V通道进行处理。由于代码基本与前面介绍的一致，这里我们不再做过多解释了。
Img = imread('couple.tiff');
hsvImg = rgb2hsv(Img);
V=hsvImg(:,:,3);
[height,width]=size(V);
V = uint8(V*255);
NumPixel = zeros(1,256);
for i = 1:height
for j = 1: width
NumPixel(V(i,j) + 1) = NumPixel(V(i,j) + 1) + 1;
end
end
ProbPixel = zeros(1,256);
for i = 1:256
ProbPixel(i) = NumPixel(i) / (height * width * 1.0);
end
CumuPixel = cumsum(ProbPixel);
CumuPixel = uint8(255 .* CumuPixel + 0.5);
for i = 1:height
for j = 1: width
V(i,j) = CumuPixel(V(i,j));
end
end
V = im2double(V);
hsvImg(:,:,3) = V;
outputImg = hsv2rgb(hsvImg);
imshow(outputImg);
最后，来对比一下不同方法对彩色图像的处理效果。下面的左图是采用R、G、B三分量分别处理得到的结果。右图是对HSV空间下V通道处理之结果。显然，右图的效果更理想，而左图则出现了一定的色彩失真。事实上，对彩色图像进行直方图均衡是图像处理研究领域一个看似简单，但是一直有人在研究的话题。我们所说的对HSV空间中V分量进行处理的方法也是比较基本的策略。很多相关的研究文章都提出了更进一步的、适应性更强的彩色图像直方图均衡化算法。有兴趣的读者可以参阅相关文献以了解更多。
分别处理R、G、B三个分量之结果 转换到HSV空间后处理V分量
先来看一下待处理的图像效果：
下面是利用CLAHE算法处理之后得到的两个效果（后面我们还会具体介绍我们所使用的策略）：
效果图A 效果图B
对于一幅图像而言，它不同区域的对比度可能差别很大。可能有些地方很明亮，而有些地方又很暗淡。如果采用单一的直方图来对其进行调整显然并不是最好的选择。于是人们基于分块处理的思想提出了自适应的直方图均衡算法AHE。维基百科上说的也比较明白：AHE improves on this by transforming each pixel with a transformation function derived from a neighbourhood region. 但是这种方法有时候又会将一些噪声放大，这是我们所不希望看到的。于是荷兰乌得勒支大学的Zuiderveld教授又引入了CLAHE，利用一个对比度阈值来去除噪声的影响。特别地，为了提升计算速度以及去除分块处理所导致的块边缘过渡不平衡效应，他又建议采用双线性插值的方法。关于算法的介绍和描述，下面这两个资源已经讲得比较清楚。
[1] https://en.wikipedia.org/wiki/Adaptive_histogram_equalization#Contrast_Limited_AHE
[2] K. Zuiderveld: Contrast Limited Adaptive Histogram Equalization. In: P. Heckbert: Graphics Gems IV, Academic Press 1994 （http://www.docin.com/p-119164091.html）
事实上，尽管这个算法原理，然而它实现起来却仍然有很多障碍。但在此之前，笔者还需说明的是，Matlab中已经集成了实现CLAHE的函数adapthisteq()，如果你仅仅需要一个结果，其实直接使用这个函数就是最好的选择。我给出一些示例代码用以生成前面给出之效果。函数adapthisteq()只能用来处理灰度图，如果要处理彩色图像，则需要结合自己编写的代码来完成。上一篇文章介绍了对彩色图像进行直方图均衡的两种主要策略：一种是对R、G、B三个通道分别进行处理；另外一种是转换到另外一个色彩空间中再进行处理，例如HSV（转换后只需对V通道进行处理即可）。
首先，我们给出对R、G、B三个通道分别使用adapthisteq()函数进行处理的示例代码：
img = imread('space.jpg');
rimg = img(:,:,1);
gimg = img(:,:,2);
bimg = img(:,:,3);
resultr = adapthisteq(rimg);
resultg = adapthisteq(gimg);
resultb = adapthisteq(bimg);
result = cat(3, resultr, resultg, resultb);
imshow(result);
下面程序将原图像的色彩空间转换到LAB空间之后再对L通道进行处理。
clear;
img = imread('space.jpg');
cform2lab = makecform('srgb2lab');
LAB = applycform(img, cform2lab);
L = LAB(:,:,1);
LAB(:,:,1) = adapthisteq(L);
cform2srgb = makecform('lab2srgb');
J = applycform(LAB, cform2srgb);
imshow(J);
上述程序所得之结果如图B所示。
如果你希望把这个算法进一步提升和推广，利用用于图像去雾、低照度图像改善和水下图像处理，那么仅仅知其然是显然不够的，你还必须知其所以然。希望我下面一步一步实现的代码能够帮你解开这方面的困惑。鉴于前面所列之文献已经给出了比较详细的算法描述，下面将不再重复这部分内容，转而采用Matlab代码来对其中的一些细节进行演示。
首先来从灰度图的CLAHE处理开始我们的讨论。为此清理一下Matlab的环境。然后，读入一张图片（并将其转化灰度图），获取图片的长、宽、像素灰度的最大值、最小值等信息。
clc;
clear all;
Img = rgb2gray(imread('space.jpg'));
[h,w] = size(Img);
minV = double(min(min(Img)));
maxV = double(max(max(Img)));
imshow(Img);
图像的初始状态显示如下。此外该图的 Height = 395，Width = 590，灰度最大值为255，最小值为8。
我们希望把原图像水平方向分成8份，把垂直方向分成4份，即原图将被划分成4 × 8 = 32个SubImage。然后可以算得每个块（tile）的height = 99，width = 74。注意，由于原图的长、宽不太可能刚好可被整除，所以我在这里的处理方式是建立一个稍微大一点的图像，它的宽和长都被补上了deltax和deltay，以保证长、宽都能被整除。
NrX = 8;
NrY = 4;
HSize = ceil(h/NrY);
WSize = ceil(w/NrX);
deltay = NrY*HSize - h;
deltax = NrX*WSize - w;
tmpImg = zeros(h+deltay,w+deltax);
tmpImg(1:h,1:w) = Img;
对长和宽进行填补之后，对新图像的一些必要信息进行更新。
new_w = w + deltax;
new_h = h + deltay;
NrPixels = WSize * WSize;然后指定图像中直方图横坐标上取值的计数（也就指定了统计直方图上横轴数值的间隔或计数的精度），对于色彩比较丰富的图像，我们一般都要求这个值应该大于128。
% NrBins - Number of greybins for histogram ("dynamic range")
NrBins = 256;
然后用原图的灰度取值范围重新映射了一张Look-Up Table（当然你也可以直接使用0~255这个范围，这取决你后续建立直方图的具体方法），并以此为基础为每个图像块（tile）建立直方图。
LUT = zeros(maxV+1,1);
for i=minV:maxV
LUT(i+1) = fix(i - minV);%i+1
end
Bin = zeros(new_h, new_w);
for m = 1 : new_h
for n = 1 : new_w
Bin(m,n) = 1 + LUT(tmpImg(m,n) + 1);
end
end
Hist = zeros(NrY, NrX, 256);
for i=1:NrY
for j=1:NrX
tmp = uint8(Bin(1+(i-1)*HSize:i*HSize, 1+(j-1)*WSize:j*WSize));
%tmp = tmpImg(1+(i-1)*HSize:i*HSize,1+(j-1)*WSize:j*WSize);
[Hist(i, j, :), x] = imhist(tmp, 256);
end
end
Hist = circshift(Hist,[0, 0, -1]);注意：按通常的理解，上面这一步我们应该建立的直方图（集合）应该是一个4×8=32个长度为256的向量（你当然也可以这么做）。但由于涉及到后续的一些处理方式，我这里是生成了一个长度为256的4×8矩阵。Index = 1的矩阵其实相当于是整张图像各个tile上灰度值=0的像素个数计数。例如，我们所得的Hist(:, :, 18)如下。这就表明图像中最左上角的那个tile里面灰度值=17的像素有零个。同理，它右边的一个tile则有46个灰度值=17的像素。
Hist(:,:,18) =
0 46 218 50 14 55 15 7
0 0 21 18 114 15 74 73
0 1 0 0 2 67 124 82
0 0 0 0 0 1 9 2
然后来对直方图进行裁剪。Matlab中内置的函数adapthisteq()中cliplimit参数的取值范围是0~1。这里我们所写的方法则要求该值>1。当然这完全取决你算法实现的策略，它们本质上并没有差异。然后我们将得到新的（裁剪后的）映射直方图。
ClipLimit = 2.5;
ClipLimit = max(1,ClipLimit * HSize * WSize/NrBins);
Hist = clipHistogram(Hist,NrBins,ClipLimit,NrY,NrX);
Map=mapHistogram(Hist, minV, maxV, NrBins, NrPixels, NrY, NrX);
因为这里没有具体给出clipHistogram函数的实现，所以此处我希望插入一部分内容来解释一下我的实现策略（也就是说，在实际程序中并不需要包含这部分）。我们以图像最左上角的一个tile为例，它的原直方图分布可以用下面代码来绘出：
tmp_hist = reshape(Hist(1,1,:), 1, 256);
plot(tmp_hist)输出结果下图中的左图所示。
如果我们给ClipLimit赋初值为2.5，则经过语句ClipLimit = max(1,ClipLimit * HSize * WSize/NrBins);计算之后，ClipLimit将变成71.54。然后我们再用上述代码绘制新的直方图，其结果将如上图中的右图所示。显然，图中大于71.54的部分被裁剪掉了，然后又平均分配给整张直方图，所以你会发现整张图都被提升了。这就是我们这里进行直方图裁剪所使用的策略。但是再次强调，matlab中的内置函数adapthisteq()仅仅是将这个参数进行了归一化，这与我们所使用的方法并没有本质上的区别。
继续回到程序实现上的讨论。最后，也是最关键的步骤，我们需要对结果进程插值处理。这也是Zuiderveld设计的算法中最复杂的部分。
yI = 1;
for i = 1:NrY+1
if i == 1
subY = floor(HSize/2);
yU = 1;
yB = 1;
elseif i == NrY+1
subY = floor(HSize/2);
yU = NrY;
yB = NrY;
else
subY = HSize;
yU = i - 1;
yB = i;
end
xI = 1;
for j = 1:NrX+1
if j == 1
subX = floor(WSize/2);
xL = 1;
xR = 1;
elseif j == NrX+1
subX = floor(WSize/2);
xL = NrX;
xR = NrX;
else
subX = WSize;
xL = j - 1;
xR = j;
end
UL = Map(yU,xL,:);
UR = Map(yU,xR,:);
BL = Map(yB,xL,:);
BR = Map(yB,xR,:);
subImage = Bin(yI:yI+subY-1,xI:xI+subX-1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sImage = zeros(size(subImage));
num = subY * subX;
for i = 0:subY - 1
inverseI = subY - i;
for j = 0:subX - 1
inverseJ = subX - j;
val = subImage(i+1,j+1);
sImage(i+1, j+1) = (inverseI*(inverseJ*UL(val) + j*UR(val)) ...
+ i*(inverseJ*BL(val) + j*BR(val)))/num;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
output(yI:yI+subY-1, xI:xI+subX-1) = sImage;
xI = xI + subX;
end
yI = yI + subY;
end这个地方，作者原文中已经讲得比较清楚了，我感觉我也没有必要狗尾续貂，班门弄斧了。下面截作者原文中的一段描述，足以说明问题。
最后来看看我们处理的效果如何（当然，这里还需要把之前我们填补的部分裁掉）。
output = output(1:h, 1:w);
figure, imshow(output, []);来看看结果吧~可以对比一下之前的灰度图，不难发现，图像质量已有大幅改善。
2833
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
