图像对比度增强算法

一、全局对比度增强$^1$

1. 直方图均衡化 Histogram Equalization

a. 基本思想

对于图像中的灰度点做映射，使得整体图像的灰度大致符合均匀分布。增强图像的对比度。
- 对于灰度值连续的情况，使用灰度的累积分布函数CDF做转换函数，可以使得输出图像的灰度符合均匀分布。
- 对于灰度值不连续的情况，存在舍入误差，得到的灰度分布大致符合均匀分布。
- 直观地理解，如果某一个灰度范围（如200-201）的像素点很少，那么它的概率密度值就会很小，所以CDF在200-201附近的增长变化就会很小；反之，如果某一个灰度范围（如100-101）的像素点很多，CDF在100-101附近的增长变化会很大。总体来看，以灰度为横轴，CDF为纵轴画曲线。这种向上凸的曲线，很像gamma变换: $s=cr^\gamma$中$\gamma<1$的情形。将灰度集中的部分拉伸，而将灰度不集中的部分压缩，达到提高对比度的效果。
- 直方图均衡可以看做自适应的gamma变换或者分段变换。前者的优势在于，不需要指定任何参数，所有运算都是基于图像本身的。

b. 算法

1. 根据图像灰度计算灰度概率密度函数$PDF$
2. 计算累积概率分布函数$CDF$
3. 将$CDF$归一化到原图灰度取值范围，如[0,255]。
4. 之后$CDF$四舍五入取整，得到灰度转换函数$s_k=T(r_k)$
5. 将$CDF$作为转换函数，将灰度为$r_k$的点转换为$s_k$灰度

c. matlab实验

代码$^2$：

%% 直方图均衡
clear all;clc;close all;
ImgFile='E:\图像处理\冈萨雷斯图片库\DIP3E_Original_Images_CH03\Fig0310(b)(washed_out_pollen_image).tif';
ImgIn=imread(ImgFile);
ImgHistEq=histeq(ImgIn,256);
figure;subplot(121);imshow(uint8(ImgIn));title('原图');
subplot(122);imshow(ImgHistEq);title('全局灰度增强 - 直方图均衡');
figure;subplot(121);imhist(ImgIn,256);
axis([0 255 0 1e5]);title('原图的直方图');
subplot(122);imhist(ImgHistEq,256);axis([0 255 0 1e5]);title('直方图均衡化后的直方图');
% 自定义直方图均衡
[counts,x]=imhist(ImgIn,256);
cum_counts=cumsum(counts);
cum_counts=uint8(cum_counts/max(cum_counts)*255);% 转化函数
figure;plot(x,cum_counts);axis([0 255 0 255]);
xlabel('原图灰度');ylabel('转换后灰度');title('原图CDF转化的灰度映射函数');
ImgOut=nan(size(ImgIn));
for i=1:length(x)
    ImgOut(ImgIn==x(i))=cum_counts(i);
end
ImgOut=uint8(ImgOut);
figure;imshow(uint8(ImgOut));title('自定义直方图均衡')
figure;imhist(ImgOut,256);axis([0 255 0 1e5]);title('自定义直方图均衡的直方图')

输出：

图1. 原图和均衡化后的图。均衡后图像对比度明显增强，纹理和边缘更加清晰。

图2. 原图和均衡化后的图的直方图。原图灰度分布在很狭窄的区间约[90,140]内。均衡化后灰度较为均匀的分布在整个[0, 255]区间内。

图3. 原图CDF归一化到[0, 255]后作为灰度变换函数。为了方便显示，画成连续曲线。

图4. 自定义均衡的直方图。和Matlab自带histeq输出基本一致。

2. 直方图匹配 Histogram Matching

a. 基本思想

我的理解:

对于图像中的灰度点做映射，使得整体图像的灰度大致符合给定的分布。
想要图像的直方图Hs和期望直方图Hd相同, 那么二者的CDF也是相同的.
也就是说, 对于某一个灰度x, CDFs(x)=CDFd(x). 所以, 如果CDFs(x)=CDFd(y). 那么我们令x=y.
也就是将原图的x灰度映射为y灰度. 相当于将图像直方图的横轴做了一个点对点的缩放操作.
这个缩放操作的结果就是图像直方图向理想直方图靠近.
在缩放过程中, 如果图像的直方图某一段灰度区间被拉长了, 因为灰度都是整数, 所以直方图的bin会变得稀疏.
所以,
最终的直方图匹配的结果就是直方图的bin变得稀疏了. 有的灰度区间没有灰度点.

• 对于某些图像，灰度分布过于不均衡，例如背景很大且均匀。直接做直方图均衡会把前景“洗白”，效果很差。
• 直方图匹配需要设定一个直方图分布。不同的图像，此分布可能不同。

b. 算法

1. 根据图像计算概率密度分布$p_r(r)$；
2. 根据$p_r(r)$计算累计分布函数$s_k=T(r_k)$；
3. 根据给定的目标分布$p_z(z)$计算累计分布函数$G(z_q)$；
4. 对于每一个$k$，找到一个$q$，使得$G(z_q)$约等于$s_k$；
5. 将原图中灰度为$k$的点变为灰度$q$；

c. Matlab实验

clear all;clc;close all;
ImgFile='E:\图像处理\冈萨雷斯图片库\DIP3E_Original_Images_CH03\Fig0323(a)(mars_moon_phobos).tif';
ImgIn=imread(ImgFile);
figure;subplot(121);imshow(ImgIn);title('原图')
subplot(122);imhist(ImgIn);title('原图Hist')
% histeq
ImgHistEq=histeq(ImgIn);
figure;subplot(121);imshow(ImgHistEq);title('Hist EQ');
subplot(122);imhist(ImgHistEq,256);title('Hist EQ')
% 根据3.35(a)估计的分布
x=[0 8 16 180 200 255];
Y=[0 70000 8000 0 3139 0];
xi=0:255;
yi = interp1q(x',Y',xi');
figure;plot(xi,yi,'.');title('目标 Hist')
ImgHistMatch = histeq(ImgIn, yi);
figure;subplot(121);imshow(ImgHistMatch);title('histeq匹配')
subplot(122);imhist(ImgHistMatch,256);title('histeq匹配直方图')

输出：

图5. 原图和直方图

图6. 直方图均衡后图像及直方图

图7. 目标直方图

图8. 直方图匹配图像和直方图

自定义直方图匹配

%% 自定义匹配
clear all;clc;close all;
ImgFile='E:\图像处理\冈萨雷斯图片库\DIP3E_Original_Images_CH03\Fig0323(a)(mars_moon_phobos).tif';
ImgIn=imread(ImgFile);
x=[0 8 16 180 200 255];
Y=[0 70000 8000 0 3139 0];
xi=0:255;
yi = interp1q(x',Y',xi');
[counts, x]=imhist(ImgIn,256);
cum_counts=cumsum(counts);% 累加CDF
% 归一化到[0 255]
cum_counts=cum_counts./numel(ImgIn);
s_k=round(cum_counts*255);
cum_yi=cumsum(yi);% 目标累加CDF
g_z=round(cum_yi./sum(yi)*255); 
figure;plot(s_k,'r');hold on;plot(g_z,'g');grid on;hold off;
legend('s_k','g_z')
axis([0 255 0 255]);title('原图累计CDF 和 目标累计CDF')
map=nan(size(s_k));
for k=1:length(s_k)
    delta=abs(g_z-s_k(k));
    [~,z]=min(delta);
    map(k)=z;
end
figure;plot(0:255,map);axis([0 255 0 255]);grid on;
xlabel('输入灰度');ylabel('输出灰度');title('灰度映射曲线')

ImgOut=zeros(size(ImgIn));
for m=0:length(map)-1
    ImgOut(ImgIn==m)=map(m+1);
end
ImgOut=uint8(ImgOut);
figure;imshow(ImgOut);title('直方图匹配')
figure;imhist(ImgOut,256);title('直方图匹配');

输出：只列出转换函数

转换函数映射图

二、局部对比度增强

1. 邻域直方图均衡

将全局直方图均衡的思想应用于邻域直方图处理中。

2. 邻域直方图匹配

将全局直方图匹配的思想应用于邻域直方图处理中。

3. 邻域统计方法

a. 基本思想

对于一幅图像，灰度的均值$M_G$表示平均光照，灰度的方差（或标准差）$\sigma_G$表示对比度的强弱。
对于图像某一个像素，邻域的灰度均值$M_L$表征邻域光照强度，邻域的方差或者标准差$\sigma_L$表示邻域对比度的强弱。
根据$M_L$和$\sigma_L$的大小自主选择想要增强的点做增强（例如乘以某一个常数$E$）。

b. 算法

1. 初始化：增强常数$E$，灰度下阈值$k_0$，标准差下阈值 $k_1$，标准差上阈值 $k_2$，窗口半宽$s$；
2. 计算图像灰度均值$M_G$和灰度标准差$\sigma_G$；
   • 对于每一个像素，计算邻域（大小为$2*step+1$的方块）内灰度均值$M_L$和标准差$\sigma_L$；
   • 如果$M_L<=k_0*M_G$并且$k_1*\sigma_G<=\sigma_L<=k_2*\sigma_G$，将像素灰度乘以 E <script type="math/tex" id="MathJax-Element-43">E</script>。

   • c. matlab代码

     测试脚本LocalEnhancement_test.m：

     %% 局部对比度增强
     clear all;clc;close all;
     ImgFile='E:\图像处理\冈萨雷斯图片库\DIP3E_Original_Images_CH03\tungsten_original.tif';
     step=1;
     para.E=4.0;
     para.k0=0.4; % 均值下阈值
     para.k1=0.02; % 标准差下阈值
     para.k2=0.4; % 标准差上阈值
     ImgIn=imread(ImgFile);
     ImgHistEq=histeq(ImgIn,256);
     ImgIn=double(ImgIn);
     % ImgIn=double(rgb2gray(ImgIn));
     [ MeanLocal,VarLocal ] = LocalStatistics( ImgIn, step );
     [ ImgOut ] = LocalEnhancement( ImgIn, MeanLocal, VarLocal,para );
     
     figure;imshow(uint8(ImgIn));title('原图');
     figure;imshow(uint8(ImgOut));title('局部统计增强');
     figure;imshow(ImgHistEq);title('全局灰度增强 - 直方图均衡');

     函数LocalEnhancement.m

     function [ ImgOut ] = LocalEnhancement( ImgIn, MeanLocal, VarLocal,para )
     %LOCALENHANCEMENT Summary of this function goes here
     % 使用局部均值和局部标准差实现局部图像增强
     %   Detailed explanation goes here
     [rows, cols]=size(ImgIn);
     MeanGlobal=mean(ImgIn(:));
     VarGlobal=std(ImgIn(:));
     for i=1:rows
         for j=1:cols
             if MeanLocal(i,j)<=para.k0*MeanGlobal...
                     &&VarLocal(i,j)>=para.k1*VarGlobal...
                     &&VarLocal(i,j)<=para.k2*VarGlobal
                 ImgIn(i,j)=para.E*ImgIn(i,j);
             end
         end
     end
     ImgOut=ImgIn;
     end
     

     函数LocalStatistics.m

     function [ MeanLocal,VarLocal ] = LocalStatistics( ImgIn, step )
     %LOCALSTATISTICS Summary of this function goes here
     % 求局部的灰度均值灰度方差
     % ImgIn - 输入灰度图，double
     % step  - 窗口大小为2*step+1，step越大，细节越少
     %   Detailed explanation goes here
     % 矩阵边缘补齐
     ImgIn=padarray(ImgIn,[step, step],'replicate');
     [rows, cols] = size(ImgIn);
     MeanLocal=zeros(size(ImgIn));
     VarLocal=zeros(size(ImgIn));
     for i = step+1:1:rows-step
         for j = step+1:1:cols-step
             LocalNeighbor=ImgIn(i-step:i+step,j-step:j+step);
             MeanLocal(i,j)=mean(LocalNeighbor(:));
             VarLocal(i,j)=std(LocalNeighbor(:));
         end
     end
     % 切割成和原图相同的大小
     MeanLocal=MeanLocal(2:end-1,2:end-1);
     VarLocal=VarLocal(2:end-1,2:end-1);
     
     end

     结果输出：
     目的是想要增强右下方的纹理，而不改变中间较亮的纹理。全局灰度增强不但不能增强右下方的纹理，而且改变了中间较亮的纹理。局部统计增强达到了目的，但是产生了一些不想要的假边缘（如左上方）和奇异点（一些散乱的白点）。
     
     图9. 原图和全局灰度增强的图。全局histeq对右下角纹理的增强效果不明显，而且改变了中间较亮区域的灰度。

     
     图10. 局部统计增强。右下方纹理得到明显的增强，而且不改变中间较亮区域的纹理。但是产生了一些边缘失真。

     ---

     2017.4.24补充:

     
     直方图均衡和其他的均衡方法都会把图像的直方图尽量调整的分布均匀一些. 尽量不会出现较高的尖峰和较大的灰度空档.这是因为这些算法不知道要增强哪些部分的灰度, 也不知道要去抑制哪些部分的灰度. 而某一部分灰度的增强必然意味着另外一部分灰度的压缩. 所以这些算法只能”雨露均沾”, 把hist中所有的灰度都调整到较为均匀的分布. 雨露均沾的结果, 可能是细节的增强和提升. 也有可能是背景被增强, 造成背景噪声放大, 还不如不增强. 针对后一种情况, 就有了CLAHE等限制对比度的算法.
     

     ---

     1. 本文主要参考冈萨雷斯《数字图像处理》第三版英文版
     2. 如果matlab读入*.tif图片出错，参见：
        http://jingyan.baidu.com/article/aa6a2c14f9441b0d4c19c4b0.html