图像处理：课程学习系列（3）——图像增强之空间域滤波（数学原理+matlab实现）

1、卷积运算定义

图像是二维且离散的。二维离散卷积运算定义为：
$$y(j,i)=h(j,i)\ast x(j,i)\to y(j,i)={\Sigma }_m{\Sigma }_{n}h(m,n)x(j+m,i+n)$$
其中，$h(m,n)$为（图像增强）系统冲击响应，也就是我们常说的卷积核矩阵，

2、低通滤波器

2.1、作用

（1）把细小的噪声去除，提取感兴趣物体

（2）图像平滑，比如去除皱纹

缺点：会让图像的边缘变模糊。

2.2、均值滤波器

均值滤波器的卷积核如下：
$$h(m,n)=\left[\begin{array}{ccc}1& 1& 1\\ 1& 1& 1\\ 1& 1& 1\end{array}\right]$$
上述的二维离散卷积运算变形为：
$$\begin{gathered} y(j,i)=\frac{1}{9}{\Sigma }_{m=-1}^1{\Sigma }_{n=-1}^{1}x(j+m.i+n) \\ \downarrow \\ y(j,i)=\frac{1}{9}[x(j-1,i-1)+x(j-1,i)+x(j-1,i+1)+x(j,i)+\cdots ] \end{gathered}$$
tips：就是一个数（如下图红点所示）用自身周边九宫格（如下图九宫格）的平均数取代；其中9=1+1+1+1+1+1+1+1+1，为矩阵的总权数。

2.3、高斯滤波器

高斯滤波器的卷积核如下：
$$h(m,n)=\left[\begin{array}{ccc}1& 2& 1\\ 2& 4& 2\\ 1& 2& 1\end{array}\right]$$
上述的二维离散卷积运算变形为：
$$y(j,i)=\frac{1}{16}[x(j-1,i-1)+2\ast x(j-1,i)+x(j-1,i+1)+4\ast x(j,i)+\cdots ]$$
tips：就是一个数用自身周边九宫格的加权平均数取代；其中16=1+2+1+2+4+2+1+2+1，为矩阵的总权数。

2.4、matlab实现

close all;
 
 A = imread('ROI-first\1_1.jpg');
 I = rgb2gray(A);
 subplot(1,3,1); 
 imshow(I);  
 title('原图'); 
 
 %均值滤波器
 h1=[1 1 1;1 1 1 ;1 1 1];
 I1 = lfilter(I,h1);
 subplot(1,3,2); 
 imshow(I1);  
 title('均值滤波器平滑后的图像'); 
 
  %高斯滤波器
 h2=[1 2 1;2 4 2 ;1 2 1];
 I2 = lfilter(I,h2);
 subplot(1,3,3); 
 imshow(I2);  
 title('高斯滤波器平滑后的图像'); 
 
 function I = lfilter(I,H)
   I = double(I);
   h_x = size(H,1);
   h_y = size(H,2);
   count = 0;
   for i = 1:h_x
     for j = 1:h_y
         count = count + H(i,j);
     end
   end
   
   x = size(I,1);
   y = size(I,2);
   new = zeros(x+2,y+2);
   new(2:end-1,2:end-1) = I;
   for i = 1:length(I(:,1))
       for j = 1:length(I(1,:))
           I(i,j) = 1/count *(H(1,1)*new(i,j) + H(1,2)*new(i,j+1) + H(1,3)*new(i,j+2) ...
               + H(2,1)*new(i+1,j) + H(2,2)*new(i+1,j+1) + H(2,3)*new(i+1,j+2)  ...
               + H(3,1)*new(i+2,j) + H(3,2)*new(i+2,j+1) + H(3,3)*new(i+2,j+2));
       end
   end   
   I = mat2gray(I);
 end

3、中值滤波器

3.1、存在的理由

低通滤波器虽然会抑制噪声，但会造成边缘模糊。

3.2、作用

（1）突出亮（暗）点。
（2）消除孤立的亮（暗）点。
（3）更接近它周边的点。
（4）能够有效地去除脉冲噪声。
（5）去除噪声的同时，比较好地保留了边缘。

3.3、工作原理

一般也是利用该像素附近的九宫格（也就是3*3，不排除用更小/大的卷积核），
（1）对九宫格的像素值排序；
（2）取出中间那个值，若是偶数，则取中间两个数的平均值；
（3）用中值取代该像素值。

3.4、matlab实现

close all;
 
 A = imread('ROI-first\1_1.jpg');
 I = rgb2gray(A);
 
 subplot(1,2,1); 
 imshow(I);  
 title('原图'); 
 
 x = size(I,1);
 y = size(I,2);
 new = zeros(x+2,y+2);
 new(2:end-1,2:end-1) = I;
 
   for i = 1:length(I(:,1))
       for j = 1:length(I(1,:))
           I(i,j) = median([new(i,j) new(i,j+1) new(i,j+2) new(i+1,j) new(i+1,j+1) new(i+1,j+2)  ...
               new(i+2,j) new(i+2,j+1) new(i+2,j+2)]);
       end
   end
   
 subplot(1,2,2); 
 imshow(I);  
 title('中值滤波器平滑后的图像'); 

4、高通滤波器

4.1、作用

（1）增强高频。
（2）突出边缘。
（3）图像锐化。
tips：主要是通过检测出图像边缘，然后与原图像叠加，进行图像增强。

4.2、一阶差分算子

边缘定义为：
$$|▽f|=|G_x|+|G_y|$$
（1）Robert算子
$$\begin{gathered} G_x=f(x+1,y+1)-f(x,y) \\ {G}_y=f(x,y+1)-f(x+1,y) \end{gathered}$$
水平和垂直方向的卷积核分别为：
$$\left[\begin{array}{cc}-1& 0\\ 0& 1\end{array}\right]\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right]$$
（2）Prewitt算子
$$\begin{gathered} G_x=(f(x+1,y+1)+f(x+1,y)+f(x+1,y-1))-((f(x-1,y+1)+f(x-1,y)+f(x-1,y-1)) \\ {G}_y=(f(x-1,y+1)+f(x,y+1)+f(x+1,y+1))-((f(x-1,y-1)+f(x,y-1)+f(x+1,y-1)) \end{gathered}$$
水平和垂直方向的卷积核分别为：
$$\left[\begin{array}{ccc}-1& 0& 1\\ -1& 0& 1\\ -1& 0& 1\end{array}\right]\left[\begin{array}{ccc}-1& -1& -1\\ 0& 0& 0\\ 1& 1& 1\end{array}\right]$$
（3）Sobel算子
$$\begin{gathered} G_x=(f(x+1,y-1)+2f(x+1,y)+f(x+1,y+1))-((f(x-1,y+1)+2f(x-1,y)+f(x-1,y-1)) \\ {G}_y=(f(x-1,y-1)+2f(x,y+1)+f(x+1,y+1))-((f(x-1,y-1)+2f(x,y-1)+f(x+1,y+1)) \end{gathered}$$
水平和垂直方向的卷积核分别为：
$$\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right]\left[\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right]$$

4.3、二阶差分算子

边缘定义为：
$$|{▽}^{2}f|=\frac{{\partial }^{2}f}{{\partial }^{2}x}+\frac{{\partial }^{2}f}{{\partial }^{2}y}$$
（1）Laplace算子
$${▽}^{2}f=f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)-4f(x,y)$$
卷积核为：
$$\left[\begin{array}{ccc}0& -1& 0\\ -1& 4& -1\\ 0& -1& 0\end{array}\right]$$

4.4、一、二阶差分算子的差异

（1）一阶导数可以检测图像中的某像素点是否在边缘上。一阶算子对噪声敏感度较低，可以检测边缘的方向，但不能判断像素点在边缘的内部还是外部。

（2）二阶导数可以判断一个边缘像素点在亮或暗的一边。二阶算子对噪声敏感（双倍加强），不能检测边缘方向，但能进行边缘定位，具有旋转不变性。

4.5、matlab实现

 close all;
 
 A = imread('ROI-first\1_1.jpg');
 I = rgb2gray(A);
 
 subplot(2,3,1); 
 imshow(I);  
 title('原图'); 
 
 threshold = 0.12; %阈值，用于过滤掉不必要的边缘
 
 %% 一阶算子
 %Sobel算子
 hx=[-1 -2 -1;0 0 0 ;1 2 1];%生产sobel垂直梯度模板
 hy=hx'; %生产sobel水平梯度模板
%方法一：
%  gradx=filter2(hx,I,'same');
%  gradx=abs(gradx); %计算图像的sobel垂直梯度
%  grady=filter2(hy,I,'same');
%  grady=abs(grady); %计算图像的sobel水平梯度
%方法二：
% grad = edge(I,‘sobel’);
 gradx = oper(I,hx)*255;
 grady = oper(I,hy)*255;
 grad=gradx+grady;  %得到图像的sobel梯度
 for i = 1:length(grad(:,1))
    for j = 1:length(grad(1,:))
        if grad(i,j) < 255*threshold %256为灰度级别数
            grad(i,j) = 0;
        end
    end
 end
 subplot(2,3,2);
 imshow(grad,[]);
 title('Sobel算子提取的图像边缘');
 
 %Robert算子
 hx = [-1 0;0 1];
 hy = [0 -1;1 0];
%  gradx=filter2(hx,I,'same');
%  gradx=abs(gradx); %计算图像的Robert垂直梯度
%  grady=filter2(hy,I,'same');
%  grady=abs(grady); %计算图像的Robert水平梯度
 gradx = oper(I,hx)*255;
 grady = oper(I,hy)*255;
 grad=gradx+grady;  %得到图像的Robert梯度
 for i = 1:length(grad(:,1))
    for j = 1:length(grad(1,:))
        if grad(i,j) < 255*threshold %256为灰度级别数
            grad(i,j) = 0;
        end
    end
 end
 subplot(2,3,3);
 imshow(grad,[]);
 title('Robert算子提取的图像边缘');
 
 %Prewitt算子
 hx=[-1 -1 -1;0 0 0 ;1 1 1];%生产Prewitt垂直梯度模板
 hy=hx'; %生产Prewitt水平梯度模板  
 gradx=filter2(hx,I,'same');
 gradx=abs(gradx); %计算图像的Prewitt垂直梯度
 grady=filter2(hy,I,'same');
 grady=abs(grady); %计算图像的Prewitt水平梯度
%  gradx = oper(I,hx)*255;
%  grady = oper(I,hy)*255;
 grad=gradx+grady;  %得到图像的Prewitt梯度
 for i = 1:length(grad(:,1))
    for j = 1:length(grad(1,:))
        if grad(i,j) < 255*threshold %255为最大灰度值
            grad(i,j) = 0;
        end
    end
 end
 subplot(2,3,4);
 imshow(grad,[]);
 title('Prewitt算子提取的图像边缘');
 
 %% 二阶算子
 % Laplace算子
%  h = [0 1 0;1 -4 1 ;0 1 0];
 h = [1 1 1;1 -8 1 ;1 1 1]; %其他常用模板1
 grad = oper(I,h)*255;
%  for i = 1:length(grad(:,1))
%     for j = 1:length(grad(1,:))
%         if grad(i,j) < 255*threshold %255为最大灰度值
%             grad(i,j) = 0;
%         end
%     end
%  end
 subplot(2,3,5);
 imshow(grad,[]);
 title('Laplace算子提取的图像边缘');
  
 %% Canny算子
 grad = edge(I,'canny');
 subplot(2,3,6);
 imshow(grad,[]);
 title('Canny算子提取的图像边缘');
 
 %% 函数区
 function I = oper(I,H)
   I = double(I);
   s = size(H,1);
   x = size(I,1);
   y = size(I,2);
   switch s
       case 2 %核卷积矩阵大小为2*2
           new = zeros(x+1,y+1);
           new(1:end-1,1:end-1) = I;
           for i = 1:length(I(:,1))
               for j = 1:length(I(1,:))
                   I(i,j) = abs(H(1,1)*new(i,j) + H(1,2)*new(i,j+1) + H(2,1)*new(i+1,j) + H(2,2)*new(i+1,j+1));
               end
           end
        case 3 %核卷积矩阵大小为3*3
           new = zeros(x+2,y+2);
           new(2:end-1,2:end-1) = I;
           for i = 1:length(I(:,1))
               for j = 1:length(I(1,:))
                   I(i,j) = abs(H(1,1)*new(i,j) + H(1,2)*new(i,j+1) + H(1,3)*new(i,j+2) ...
                       + H(2,1)*new(i+1,j) + H(2,2)*new(i+1,j+1) + H(2,3)*new(i+1,j+2)  ...
                       + H(3,1)*new(i+2,j) + H(3,2)*new(i+2,j+1) + H(3,3)*new(i+2,j+2));
               end
           end   
   end
   I = mat2gray(I);
 end

tips：上述代码中包含：
（1）三种方法（edge函数，卷积核函数，自定义的函数）求解边缘算子；
（2）阈值过滤；
（3）区间转换。
（4）还增加了canny算子