基于直方图的图像增强算法（HE、CLAHE、Retinex）之（二）

信息。

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);