matlabahe，基于直⽅图的图像增强算法（HE、CLAHE、
Retinex）之（⼆）
作为图像增强算法系列的第⼆篇⽂章，下⾯我们将要介绍功能强⼤、⽤途⼴泛、影响深远的对⽐度有限的⾃适应直⽅图均衡
(CLAHE，Contrast Limited Adaptive Histogram Equalization)算法。尽管最初它仅仅是被当作⼀种图像增强算法被提出，但是现今在图像去雾、低照度图像增强，⽔下图像效果调节、以及数码照⽚改善等⽅⾯都有应⽤。这个算法的算法原理看似简单，但是实现起来却并不那么容易。我们将结合相应的Matlab代码来对其进⾏解释。希望你在阅读本⽂之前对朴素的直⽅图均衡算法有所了解，相关内容可以参见本系列的前⼀篇⽂章：
先来看⼀下待处理的图像效果：
下⾯是利⽤CLAHE算法处理之后得到的两个效果(后⾯我们还会具体介绍我们所使⽤的策略)：
效果图A                                                                        效果图B
对于⼀幅图像⽽⾔，它不同区域的对⽐度可能差别很⼤。可能有些地⽅很明亮，⽽有些地⽅⼜很暗淡。如果采⽤单⼀的直⽅图来对其进⾏调整显然并不是最好的选择。于是⼈们基于分块处理的思想提出了⾃适应的直⽅图均衡算法AHE。上说的也⽐较明⽩：AHE improves on this by transforming each pixel with a transformation function derived from a neighbourhood region. 但是这种⽅法有时候⼜会将⼀些噪声放⼤，这是我们所不希望看到的。于是荷兰乌得勒⽀⼤学的Zuiderveld教授⼜引⼊了CLAHE，利⽤⼀个对⽐度阈值来去除噪声的影响。特别地，为了提升计算速度以及去除分块处理所导致的块边缘过渡不平衡效应，他⼜建议采⽤双线性插值的⽅法。关于算法的介绍和描述，下⾯这两个资源已经讲得⽐较清楚。
事实上，尽管这个算法原理，然⽽它实现起来却仍然有很多障碍。但在此之前，笔者还需说明的是，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));
matlab直方图%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);
最后，也是最关键的步骤，我们需要对结果进程插值处理。这也是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);

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。