图像处理：课程学习系列（2）——图像增强之灰度变换+代数运算（数学原理+matlab实现）

1、为什么要图像增强

（1）视觉效果不佳
（2）噪声污染
（3）难以分析理解

2、定义

按照特定的需要突出或去除图像中的某些信息。

3、图像增强的方法

（1）灰度变化

分为线性变换和非线性变换。

（1.1）线性变换

$s=T[r]$：有反转操作，线性函数加强，以及分段函数变换，对应的数学公式分别是$$\begin{gathered} s=255-r \\ s=Ar+B \\ s=\{\begin{array}{rcl}{T}_1[r]& & 0<r<a\\ T_2[r]& & a\le r<b\\ T_3[r]& & b\le r<255\end{array} \end{gathered}$$

（1.2）非线性变换

包含对数变换、幂次变换和直方图均衡化等。
（1.2.1）对数变换：对比较暗的区域进行拉伸变量，而对较亮的区域则没有太多的处理，比较适合窄带低灰度的图像。
$$s=c\ast \log (1+r)$$
（1.2.2）幂次变换：可以根据图像的特性选取不同的$\gamma$从而达到不同的效果，当$\gamma >1$时，更多的是针对高灰度进行拉伸，使图像变暗；当$\gamma <1$时，更多的是针对低灰度进行拉伸，使图像变亮。
$$s=c\ast r^{\gamma }$$


（1.2.2）直方图均衡化：
前提：灰度级$D_A$，$D_B$；变换函数$D_B=f(D_A)$；灰度直方图$H(D_A)$，$H(D_B)$。

由上图阴影面积相等可得，
$${\int }_{D_B}^{D_B+\Delta D_B}{H}_B(D_B)dD={\int }_{D_A}^{D_A+\Delta D_A}{H}_A(D_A)dD$$
当$\Delta D_B$和$\Delta D_A$比较小的时候，上式可写为，
$$H_B\Delta D_B=H_A\Delta D_A$$
当$\Delta D_B$和$\Delta D_A$趋于0时，又可得，
$$H_B(D_B)=\frac{H_A(D_A)}{\frac{d{D}_B}{d{D}_A}}=\frac{H_A(D_A)}{f^{\prime }(D_A)}$$
求得的$H_B(D_B)$预期为一个常数，有
$$f^{\prime }(D_A)=\frac{D_m}{A_0}\cdot H_A(D_A)$$
两边积分有，
$$D_B=f(D_A)=\frac{D_m}{A_0}\cdot \sum _0^{D_A}{H}_A(D_A)$$
其中，$D_m$为灰度级，$A_0$为图像像素总数。根据上式就能求得均衡化后的图像灰度直方图。

（2）代数运算

（2.1）加法运算

作用：去除叠加性噪声。
公式：$C(x，y)=A(x,y)+B(x,y)$
描述：通过对同一幅图像多次叠加求取平均值。

（2.2）减法运算

作用：分割特定区域；检测场景变化。
公式：$C(x，y)=A(x,y)-B(x,y)$

（2.3）乘法运算

作用：获取图像中的特定区域。
公式：$C(x，y)=A(x,y)\ast B(x,y)$

（2.4）除法运算

作用：校正成像设备的非线性影响；检测两幅图像间的区别。
公式：$C(x，y)=A(x,y)/B(x,y)$

（3）空间域滤波

（4）频域滤波

tips：
【1】（1）-（3）直接对图像中像素的灰度级进行操作，属于空间域增强，（4）对图像先进行傅里叶变化转化成频域，再通过滤波器过滤一定频率的信号。
【2】 灰度归一化的作用：把图像的灰度区间大小限制为[0，255]，不用担心会超标。
【3】 直方图均衡化的代码思路：
根据公式和累计灰度计算变化前后的灰度级对应关系（这里可以用两个0-255的一维向量建立对应关系），转换后即可得新图，也可以用matlab自带的直方图均衡化代码（help histeq可查看详细情况）。

4、matlab实现

 close all
%  clear all
 
 A = imread('1.jpg');
 
 g = rgb2gray(A);
 g1 = 255 - g;
 g2 = 1.2*g + 1;
 
 %对比度拉伸
 a1 = 0.8;
 b1 = 0;
 a2 = 1.2;
 b2 = 0;
 a3 = 0.8;
 b3 = 0;
 g3 = g;
 for i = 1:size(g3,1)
     for j = 1:size(g3,2)
         if g3(i,j) < round(255/3)
            g3(i,j) = a1* g3(i,j) + b1;
         elseif g3(i,j) >= round(255/3) && g3(i,j) < round(255/3*2)
             g3(i,j) = a2* g3(i,j) + b2;
         else
             g3(i,j) = a3* g3(i,j) + b3;
         end
     end
 end
 
 g4 = mat2gray(2*log(1+double(g)));
 g5 = mat2gray(g);
 g6 = mat2gray(double(g).^2);
 g7 = myHistogramEqualization(g);
 
 figure;
 subplot(2,4,1);
 imshow(g),title('原始图像');
 subplot(2,4,2);
 imshow(g1),title('反转y=255-y');
 subplot(2,4,3);
 imshow(g2),title('线性变换y=1.2x+1');
 subplot(2,4,4);
 imshow(g3),title('分段函数变换（对比度拉伸）');
 subplot(2,4,5);
 imshow(g4),title('对数变换y=2*log（1+y）');
 subplot(2,4,6);
 imshow(g5),title('灰度归一化');
 subplot(2,4,7);
 imshow(g6),title('幂次变换');
 subplot(2,4,8);
 imshow(g7),title('直方图均衡化');
 
 function img2 = myHistogramEqualization(img)
    img1 = double(img);
    [r,c] = size(img1);%获取图像的高r和宽c
    %统计图像中每个灰度级出现的次数
    count = zeros(1,256);
    for i=1:r
        for j=1:c
            count(1,img(i,j)+1) = count(1,img(i,j)+1)+1;
        end
    end
    %统计图像中每个灰度级出现的概率
    p = zeros(1,256);
    for i=1:256
        p(1,i) = count(1,i)/(r*c);
    end
    img2 = im2uint8(ones(r,c));%创建一个r X c大小的1矩阵
    
    func_T = zeros(1,256);%变换函数
    p_sum = 0;
    %求直方图均衡化的变换函数
    for k = 1:256
        p_sum = p_sum + p(k);%求每个灰度级的概率之和
        func_T(k) = (256-1)*p_sum;%根据变换函数的公式求和
    end
    
    func_T_z =  round(func_T);%对变换函数进行取整
    %完成每个像素点的映射
    for i = 1:256
        findi = find(func_T_z==i);%找到灰度级为i的概率和
        len = length(findi);
        for j=1:len
            findj = find(img==(findi(j)-1));%进行对应每个像素点的映射
            img2(findj) = i;
        end
    end
end
 

效果图如下：