作为图画增强算法系列的第二篇文章，下面咱们即将介绍功能强大、用处广泛、影响深远的比照度有限的自适应直方图均衡(CLAHE，Contrast Limited Adaptive Histogram Equalization)算法。虽然开端它只是是被当作一种图画增强算法被提出，可是如今在图画去雾、低照度图画增强，水下图画作用调理、以及数码照片改进等方面都有运用。这个算法的算法原理看似简略，可是完结起来却并不那么简单。咱们将结合相应的Matlab代码来对其进行解说。期望你在阅览本文之前对朴素的直方图均衡算法有所了解，相关内容能够拜见本系列的前一篇文章：根据直方图的图画增强算法(HE、CLAHE、Retinex)之(一)
先来看一下待处理的图画作用：
下面是运用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);
上述程序之成果作用图A所示。
下面程序将原图画的颜色空间转化到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为例，它的原直方图散布能够用下面代码来绘出：
[plain] view plain copytmp_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, []);
来看看成果吧~能够比照一下之前的灰度图，不难发现，图画质量已有大幅改进。
