图像直方图均衡化和局部增强处理

一、基本原理

1.1直方图均衡化

（一提高图像的对比度，二使像素值几乎成均匀分布的图像0表示黑，1表示白，中心思想是把原始图像的灰度直方图从比较集中的某个灰度区间变成在全部灰度范围内的均匀分布，增加许多图像的局部对比度，亮度可以更好的在直方图上分布。）

    灰度级的直方图就是反映一幅图像中的灰度级与出现这种灰度的概率之间的关系的图形。设变量代表图像中的像素灰度级。在图像中，像素的灰度级可作归一化处理，值将限定在下述范围之内。对于一幅给定的图像来说，每一个像素取得[0, 1]区间内的灰度级是随机的，也就是说是一个随机变量。假定对每一瞬间它们是连续的随机变量，那么，就可以用概率密度函数来表示原始图像的灰度分布。如果用直角坐标系的横轴代表灰度级，用纵轴代表灰度的概率密度函数，这样就可以针对一幅图像在这个坐标系中作一曲线，这条曲线就是概率论中的分布密度曲线。
    为了有利于数字图像处理，必须引入离散形式。在离散形式下，用代表离散灰度级，用代表，并且有下式成立：
式中为图像中出现这种灰度的像素数，是图像中像素的总数，而就是概率论中所说的频数，在直角坐标系中作出与关系的图形，这个图形称为直方图。一幅给定的图像的灰度级分布在范围内。可以对[0, 1]区间内的任一个值进行如下变换：
也就是说，通过上述变换，每一个原始图像的像素灰度值都对应产生一个值。变换函数应满足下列条件：
1）在0≤r≤1区间内，为单值单调增加。
2）对于0≤r≤1，有0≤T(r)≤1。
    这里的第一个条件保证了灰度级从白到黑次序不变。第二个条件则保证了映射后的像素灰度值在允许的范围内。从到的反变换关系可表示为： 
    由概率论可知，如果已知随机变量的概率密度为，而随机变量是的函数，即，的概率密度可以由求出。
    因为是单调增加的，由数学分析可知，它的反函数也是单调函数。在这种情况下，当且仅当时发生，所以可以求得随机变量的分布函数为：
对上式两边求导，即可以得到随机变量的分布密度函数，为：
    可见，通过变换函数可以控制原图像灰度级的概率密度函数，从而改善原图像的灰度层次。直方图均衡化处理是以累积分布函数变换法为基础的直方图修正法。假定变换函数为：
式中是积分变量，而是的累积分布函数（CDF）。这里，累积分布函数是的函数，并且单调地从0增加到1，所以这个变换函数满足于在0≤r≤1区间内单值单增，对于0≤r≤1，有0≤T(r)≤1的两个条件。
    对式（6）求导，则：
    把结果代入式（7），则：
    由上面的推导可见，在变换后的变量的定义域内的概率密度是均匀分布的。由此可见，用累积分布函数作为变换函数可以产生一幅灰度级分布具有均匀概率密度的图像。其结果扩展了像素取值的动态范围。
    上面的修正方法是以连续随机变量为基础进行讨论的。为了对数字图像进行处理，必须引入离散形式的公式。当灰度级是离散值的时候，可以用频数近似代替概率值，即
式中是灰度级的总数目，是取这种灰度的概率，为图像中出现这种灰度的次数，是图像中像素的总数。通常把为得到均匀直方图的图像增强技术叫做直方图均衡化处理。式（6）的离散形式可表示为：
其反变换为：

1.2局部直方图处理

      在某种意义上，像素被基于整幅图像的灰度分布的变换函数修改。虽然这种全局方法适用于整个图像的增强，但存在这样的情况，增强图像中小区域的细节也是需要的。这些区域中，一些像素的影响在全局变换的计算中可能被忽略了，因为全局变换没有必要保证期望的局部增强。解决方法是以图像中每个像素的邻域中的灰度分布为基础设计变换函数。
这种直方图处理技术很容易适应局部增强。该过程是定义一个邻域，并把该区域的中心从一个像素移至另一个像素。在每一个位置，计算邻域中的点的直方图，并且得到的不是直方图均衡化，就是规定化变换函数。这个函数最终用于映射邻域中心像素的灰度。然后，邻域中心的被移至一个相邻像素位置，并重复该过程。当邻域进行逐像素平移时，由于只有邻域中的一行或一列改变，所以可在每一步移动中，以新数据更新前一个位置得到的直方图。这种方法与区域每移动一个像素位置就计算邻域中所有像素的直方图相比有明显的优点。

1.3使用局方图统计的局部增强

二、算法步骤

2.1直方图均衡化

Step1. 统计各灰度数目，共256个灰度级
Step2. 对应灰度值像素点数量增加一
Step3. 计算灰度分布密度
Step4. 计算累计直方图分布
Step5. 累计分布取整

Step6. 对灰度值进行映射（均衡化）

2.2 局部增强

（利用图像的局部区域，如局部均值，方差。梯度获取不同区域的差异情况，从而对图像不同区域进行不同的增值。常用直方图变化。）

三、结果分析

因为图像有相同的内容，直方图均衡导致的对比度增强足以补偿图像在视觉上难以区分灰度级的差别。在给出原始图像间的重要对比度差别后，说明了直方图均衡作为自适应对比度增强工具的强大作用。

从上图可以看到，上图中增强图像的放大图像是放大f(x,y)的值至范围之外形成g(x,y)而构成的一幅图像，相对于黑暗部分的是1，相对于明亮部分的是E，暗点表示在该坐标处增强过程没有放大像素值，亮点表示产生了一个增强的像素值。该图像是二值图像。

局部增强能够增强暗色区域，同时尽可能保留明亮区域不变。从图中可以看出，原始图像中央的钨丝及其支架很清晰但图像的右侧有另一根钨丝的结构，通过局部增强使其细节变得清楚。

四、完整程序

clear all
clc
A = imread('C:\Users\Administrator\Desktop\123.bmp');  
figure  
subplot(221);  
I=rgb2gray(A);
imshow(I);
title('原始图像');  
[height,width] = size(I);  
subplot(222); 
imhist(I);
title('原始图像直方图');  
  %进行像素灰度统计;  
NumPixel = zeros(1,256);%统计各灰度数目，共256个灰度级  
for i = 1:height  
    for j = 1: width
        NumPixel(I(i,j) + 1) = NumPixel(I(i,j) + 1) + 1;%对应灰度值像素点数量增加一  
    end  
end  
%计算灰度分布密度  
ProbPixel = zeros(1,256);  
for i = 1:256  
    ProbPixel(i) = NumPixel(i) / (height * width * 1.0);  
end  
%计算累计直方图分布  
CumuPixel = zeros(1,256);  
for i = 1:256  
    if i == 1  
        CumuPixel(i) = ProbPixel(i);  
    else  
        CumuPixel(i) = CumuPixel(i - 1) + ProbPixel(i);  
    end  
end  
%累计分布取整  
CumuPixel = uint8(255 .* CumuPixel + 0.5);  
%对灰度值进行映射（均衡化）  
for i = 1:height  
    for j = 1: width 
        I(i,j) = CumuPixel(I(i,j)+1);
    end  
end    
subplot(223);  
imshow(I);
title('直方图均衡化后的灰度图像');  
subplot(224);  
imhist(I);
title('均衡化后的直方图');
% Igray:灰度图像
% k0,k1：可调参数
% Ibw：提取得到数字区域
Igray=imread('C:\Users\Administrator\Desktop\10.bmp');
%=rgb2gray(a);
k0=0.3;
k1=0.01;
E=4;
masksize = 3; % 邻域范围，奇数
exsize = floor(masksize/2); % 原图片需要填充的区域
Iex = padarray(Igray, [exsize exsize], 'replicate', 'both'); % 图片填充
Is = zeros(size(Igray)); % 局部标准差
for i = 1:size(Igray, 1)
for j = 1:size(Igray, 2)
subdomain = Iex(i:i+2*exsize, j:j+2*exsize);
localmean = sum(subdomain(:))/masksize^2;
Is(i, j) = sqrt(sum((subdomain(:)-localmean).^2))/masksize;
end
end
% 确定不同区域的放大系数
meanIm = mean(Igray(:));
minIm = min(Igray(:));
meanIs = mean(Is(:));
maxIs = max(Is(:));
Im1 = Igray<(minIm+k0*(meanIm-minIm)); % 均值阈值提取
Is1 = Is>(meanIs + k1*(maxIs-meanIs)); % 标准差阈值提取
s = strel('disk',8); % 8需要根据实际图片大小来调
Is1 = imclose(Is1, s); % 标准差阈值提取闭运算
Ibw = Im1 & Is1; % 局部增强
figure(1)
imshow(Ibw);title('增强图像的放大图像');
Ibw2=Ibw.*(E-1)+1;
for i = 1:size(Igray, 1)
for j = 1:size(Igray, 2)
P(i,j)=Ibw2(i,j).*Igray(i,j);
end
end
figure(2)
imshow(P);title('增强的SEM图像');
figure(3)
imshow(Igray);title('局部增强原始图像');