直方图均衡化（HE）是一种很常用的直方图类方法，基本思想是通过图像的灰度分布直方图确定一条映射曲线，用来对图像进行灰度变换，以达到提高图像 对比度的目的。该映射曲线其实就是图像的累计分布直方图（CDF）（严格来说是呈正比例关系）。然而HE是对图像全局进行调整的方法，不能有效地提高局部 对比度，而且某些场合效果会非常差。如：
上述原图和HE结果图的直方图分别为：
因为从原图的直方图中求取的映射函数（CDF）形状为：
将它作用于原图像会导致直方图被整体右移，没有充分利用整个灰度动态范围。
为了提高图像的局部对比度，有人提出将图像分成若干子块，对子块进行HE处理，这便是AHE（自适应直方图均衡化），使用AHE处理上图得到：
结果直方图：
可 以看出结果图像的灰度较好地分布在了全部动态范围上。从结果图像上也可以看出，局部对比度的确得到了提高，视觉效果要优于HE。但是仍然有个问题：AHE 对局部对比度提高过大，导致图像失真。看看背景区，本来的黑色背景现在已经变成白色了，原因是因为背景区中的局部子块统计得到的直方图在0灰度处幅值太高 （实际上全黑子图基本上就集中在0灰度处），这样导致映射曲线斜率过高，将所有灰度值都映射到整个灰度轴的右侧，所以结果图中背景偏白（另外局部对比度过 高还会放大图像中的噪声，不过上图并没有体现这一点）。
为了解决这个问题，我们必须对局部对比度进行限制，这就是我们今天的主题：CLAHE。
从HE中我们知道，映射曲线T与CDF关系为：（M为最高灰度值，N为像素个数）
限制对比度，其实就是限制CDF的斜率，又因累计分布直方图CDF是灰度直方图Hist的积分：
反过来：
也就是说限制CDF的斜率就相当于限制Hist的幅度。
因此我们需要对子块中统计得到的直方图进行裁剪，使其幅值低于某个上限，当然裁剪掉的部分又不能扔掉，我们还需要将这部分裁剪值均匀地分布在整个灰度区间上，以保证直方图总面积不变，如下图：
可以看到，这时直方图又会整体上升了一个高度，貌似会超过我们设置的上限。其实在具体实现的时候有很多解决方法，你可以多重复几次裁剪过程，使得上升的部分变得微不足道，或是用另一种常用的方法：
设 裁剪值为ClipLimit，求直方图中高于该值的部分的和totalExcess，此时假设将totalExcess均分给所有灰度级， 求出这样导致的直方图整体上升的高度L=totalExcess/N，以upper= ClipLimit-L为界限对直方图进行如下处理：
(1)若幅值高于ClipLimit，直接置为ClipLimit；
(2)若幅值处于Upper和ClipLimit之间，将其填补至ClipLimit；
(3)若幅值低于Upper，直接填补L个像素点；
经过上述操作，用来填补的像素点个数通常会略小于totalExcess，也就是还有一些剩余的像素点没分出去，这个剩余来自于(1)(2)两处。这时我们可以再把这些点均匀地分给那些目前幅值仍然小于ClipLimit的灰度值。
这里给出一段代码：（摘自Matlab的adapthisteq.m），描述的就是上述过程：
% total number of pixels overflowing clip limit in each bintotalExcess = sum(max(imgHist - clipLimit,0)); % clip the histogram and redistribute the excess pixels in each binavgBinIncr = floor(totalExcess/numBins);upperLimit = clipLimit - avgBinIncr; % bins larger than this will be % set to clipLimit% this loop should speed up the operation by putting multiple pixels% into the "obvious" places firstfor k=1:numBins if imgHist(k) > clipLimit imgHist(k) = clipLimit; else if imgHist(k) > upperLimit % high bin count totalExcess = totalExcess - (clipLimit - imgHist(k)); imgHist(k) = clipLimit; else totalExcess = totalExcess - avgBinIncr; imgHist(k) = imgHist(k) + avgBinIncr; end endend % this loops redistributes the remaining pixels, one pixel at a timek = 1;while (totalExcess ~= 0) %keep increasing the step as fewer and fewer pixels remain for %the redistribution (spread them evenly) stepSize = max(floor(numBins/totalExcess),1); for m=k:stepSize:numBins if imgHist(m) < clipLimit imgHist(m) = imgHist(m)+1; totalExcess = totalExcess - 1; %reduce excess if totalExcess == 0 break; end end end k = k+1; %prevent from always placing the pixels in bin #1 if k > numBins % start over if numBins was reached k = 1; endend
CLAHE和AHE中另一个重要的问题：插值。
将图像进行分块处理，若每块中的像素点仅通过该块中的映射函数进行变换，则会导致最终图像呈块状效应：
为了解决这个问题，我们需要利用插值运算，也就是每个像素点出的值由它周围4个子块的映射函数值进行双线性插值得到，如下图：
上图中，为了求蓝色像素点处的值，需要利用它周围四个子块的映射函数分别做变换得到四个映射值，再对这四个值做双线性插值即可。
当 然对于边界处的像素点则不是通过四个子块进行插值，如上图红色像素点直接以一个子块的映射函数做变换，绿色像素则以两个子块做映射函数做线性插值。这里讲 的边界处像素是指落在图像左上角，左下角、右上角，右下角的四个子块中心像素点围成的四边形之外的像素。如下图，将图像分为8x8子块，边界像素即落在灰 色区域的像素点。
利用插值对图像进行处理的整体架构如下：
for (Y = 0; Y <= TileY; Y++) //TileY为Y方向网格数 { if (Y == 0) { SubY = TileYDim >> 1; YU = 0; YB = 0; } else if (Y == TileY) { SubY = TileYDim >> 1; YU = TileY-1; YB = YU; } else { SubY = TileYDim; YU = Y - 1; YB = Y; } for (X = 0; X <= TileX; X++) //TileX为X方向网格数 { if (X == 0) { SubX = TileXDim >> 1; XL = 0; XR = 0; } else if (X == TileX) { SubX = TileXDim >> 1; XL = TileX - 1; XR = XL; } else { SubX = TileXDim; XL = X - 1; XR = X; } MapLU = &pMapArray[numBins * (YU * TileX + XL)];//左上角映射函数 MapRU = &pMapArray[numBins * (YU * TileX + XR)];//右上角映射函数 MapLB = &pMapArray[numBins * (YB * TileX + XL)];//左下角映射函数 MapRB = &pMapArray[numBins * (YB * TileX + XR)];//右下角映射函数 Interpolate(pImPointer,Stride,Channel,MapLU,MapRU,MapLB,MapRB,SubX,SubY,aLUT);//插值 pImPointer += SubX; } pImPointer+=(SubY-1)*Stride; }
注意的是，上述循环需要（TileX+1）*（TileY+1）次，而不是TileX*TileY次。，原因很简单，以X方向为例，两侧边界处的半宽子块（灰色区）也各需要处理一次，如下图：
通过双线性插值可以基本消除块状效应：
对 于彩色图像，三通处理分开处理会导致严重的偏色，故我们可以将其进行颜色空间转换（如RGB转为HSV），然后仅对亮度分量处理，再反变换回RGB空间。 不过网上有高手将R、G、B统一处理[2]（也就相当于把一个像素点拆成三个像素点），这样得到的效果也不错，而且省去了颜色空间转换的时间，我们这里也 仿照他来吧：
另外，CLAHE对雾天图像处理效果也不错：
至于编程，我基本就是翻译adaphisteq.m，另外还有一些参考资源：
这里给出编译好的文件，有兴趣的朋友可以下载看看：
参考：
[1]
[2]