【学习笔记】Matlab自编高斯平滑器+Sobel算子求导

 这次准备先对函数进行封装，然后再写测试脚本

灰度化函数封装

  之前的笔记写了相关的算法，给出了三种灰度化的实现方案。但是一般情况下我们都不使用循环来遍历，而是使用切片的方式。所以这一次封装灰度化算法将会变成非常简单。
  在matlab中，函数的定义使用的是function。写完之后保存文件，就变成了一个.m的文件，这个就是你的函数文件，调用的时候要放在你测试脚本的同一个根目录下，否则会报错。

%灰度化函数
%平均值灰度化
function result = toGray(img)
    img = double(img);  %转换矩阵的数据类型为double
    gray_img = (img(:,:,1)+img(:,:,2)+img(:,:,3))./3;
    result = uint8(gray_img);
end

  ./在matlab中就是用于矩阵运算的，就是矩阵里面的每一个值都除以一个数，所以上面的操作就是将三个通道的二维矩阵提取出来相加之后得到一个矩阵，然后矩阵内所有值除以3。而function后面的result就是输出变量，相当于其它语言中的return result，就是传出result这个值。另外toGray就是函数的名字。
  在介绍灰度化算法的时候也说过了，RGB图像读取出来有三个尺寸数据，对应的分别是：行，列，维度。

卷积函数的封装

  还是之前的文章，将代码封装一下就变成了这样：

%卷积函数
function result = myconv(kernel,img)
    [k,num] = size(kernel);
    %判断传入的是否为double类型的数据，否则进行转化
    if ~isa(img,'double')
        p1 = double(img);
    else
        p1 = img;
    end
    [m,n] = size(p1);
    result = zeros(m-k+1,n-k+1);
    for i = ((k-1)/2+1):(m-(k-1)/2)
        for j = ((k-1)/2+1):(n-(k-1)/2)
            result(i-(k-1)/2,j-(k-1)/2) = sum(sum(kernel .* p1(i-(k-1)/2:i+(k-1)/2,j-(k-1)/2:j+(k-1)/2)));
        end
    end
end

高斯平滑器的封装

  首先来介绍一下高斯函数，其公式如下：

  这个是二维的高斯分布函数，是从一维的高斯函数推过来的，然后就没有什么然后了，反正按照公式敲代码就好，如果想要具体了解这个算法，可以参考一下别人的，讲的很详细。
  现在我们只需要知道高斯函数的作用。高斯函数，一维形式的高斯函数画出来其实就是高中学的正态分布曲线，所以有两个比较重要的参数就是：一个是方差sigma，另一个就是曲线尖峰的坐标（x轴）。
  现在转化为二维了，所以实际画出来的应该是处于三维空间，那么尖峰的坐标就变成了(x0,y0)，而方差依旧是sigma没变。
  而图像中的高斯平滑器，其实就是制造一个高斯分布的卷积核，然后用这个卷积核与原始图像进行卷积。这个操作就是所谓的高斯模糊。

%高斯平滑滤波器
function result = gaussian(img,k)
    gray_img = toGray(img); %图像灰度化
    %卷积核设定
    kernel = zeros(k,k);
    %sigma大小的确定
    %sigma = 0.8;
    sigma = 0.3*((k-1)*0.5-1)+0.8;
    center = 3;
    for i =1:k
        for j = 1:k
            %高斯函数计算模板参数
            %高斯二维分布公式
            temp = exp(-((i-center)^2+(j-center)^2)/(2*sigma^2));
            kernel(i,j) = temp/(2*pi*sigma^2);
        end
    end
    %归一化
    sums = sum(sum(kernel));
    kernel = kernel./sums;
    %
    %
    result = myconv(kernel,gray_img);
end

参数说明：其中的k就是卷积核的大小，可以是3也可以是5甚至更大，卷积核越大，模糊度越高。
当k=3时，其卷积核如下：

0.0009	0.0089	0.0195
0.0089	0.0929	0.2030
0.0195	0.2030	1.4434

然后再看看卷积核大小不同对图像的影响。

  可以看到，k=3时影响不大，但是k=11的时候就很明显了。那么高斯模糊有什么用呢？为什么要将图像模糊化？其最主要的就是为了“去噪”。因为图像中会存在一些噪声，在做边缘检测的时候，会误将这些噪声也当作是边缘来处理，这就不是我们希望见到的。

利用Sobel算子对图像进行求导

  这里可能会有疑问，为什么要对图像进行求导？而且图像实际上不是矩阵来的吗？所以数值都是离散化的啊，根本不满足求导函数的要求：只有连续的函数才能进行求导。
  求导得到的是什么？得到的其实就是函数的斜率。而在图像中，边缘的特征是什么？就是灰度值的突然变化。所以导数可以很好的描述灰度值变化的大小，因此我们可以利用图像求导来检测边缘。当然，这个方法非常粗糙。
  因为图像是离散化的一堆数值，所以没办法使用求导函数对图像进行求导，因此才会有这些“算子”的存在。Sobel算子求导的方式就是用一个3x3的卷积核来进行检测，这个卷积核分为x方向和y方向的。

  然后就是卷积操作了，先用x方向的卷积核对图像进行卷积，得到的就是竖直方向上的求导结果，y方向就是横向的。完了之后我们再进行融合操作。    当然，这么算有点慢，又是平方又是开根号的。所以我们还有一个近似值的操作，那就是求其绝对值。    接下来就是代码部分了，搞完之后再来看看效果。   这里是只利用了导数的大小做一个阈值操作，其实还可以靠梯度来检测，这样精度会高一点。这里的检测原理就是，当导数值大于一个阈值就判断是为边界。

%清理掉前面的数据
clc;
clear;
%导入图片，还是祭出我们的lena大姐
img = imread('photo/lena.png');
​
%高斯模糊
gray_img = gaossian(img,7);
%sobel算子
x = [1,2,1;0,0,0;-1,-2,-1];
y = [-1,0,1;-2,0,2;-1,0,1];
%求导计算
%分别得到竖直方向上的求导和横向的求导结果
gx = myconv(x,gray_img);
gy = myconv(y,gray_img);
G = abs(gx)+abs(gy);
[m,n] = size(G);
%定义一个模板
bw = zeros(m,n);
%设定阈值
num=50;
for i = 1:m
    for j = 1:n
        %将导数值高过某一个阈值的视为边缘
        if (G(i,j)>num)
            bw(i,j) = 255;
        else
            bw(i,j) = 0;
        end
    end
end
​
figure();
subplot(221);
imshow(uint8(gx));
subplot(222);
imshow(uint8(gy));
subplot(223);
imshow(uint8(G));
subplot(224);
imshow(uint8(bw));

lena大姐的图像看上去有些不太明显，我们换一张明显一点的，这张墙面的照片。


  可以看到，这一章图片出现了严重的噪声，就是那些白点。想要去除的话，一个是增大高斯卷积核的大小，另一个就是重新设定阈值，把阈值调高。

  这个是只把阈值调到90的结果，可以看到一些边缘丢失了，但还是存在一些噪声，这时候，可以利用其他的去噪方式去噪。比如这些白点属于椒盐噪声，可以使用中值滤波的方式去除。或者使用腐蚀膨胀套装，先腐蚀再膨胀也可以。