Python实现高斯平滑与Sobel边缘检测

时间：2021-05-24 01:53:40 阅读：0 评论：0 收藏：0 [点我收藏+]

标签：尺寸 ret add range read padding numpy cli 代码实现

平滑滤波与边缘检测是图像处理中非常基础与重要的部分。平滑滤波器主要有均值滤波，中值滤波，高斯滤波与双边滤波等，边缘检测主要有Sobel算子，Laplace算子，Canny算子等。本文主要就高斯滤波与Sobel算子进行原理上的介绍，并用Python进行实现。

第一部分，高斯滤波

原理

高斯滤波是一种线性滤波器，能够较好地平滑与抑制图像噪声，与均值中值滤波一样，高斯滤波也是对图像像素进行平均的一个过程，但不同的一点是，高斯滤波对图像像素进行的是加权平均，每一个像素点的值，都有其邻域内的其他像素点灰度值加权平均得到。在图像处理领域的高斯滤波其本质上是一种高通滤波器，过滤掉图像中的细节部分（高频部分），保留图像中的平滑部分（低频部分）。

正如上文说的，高斯滤波处理在图像上反映出来的是一个加权平均的过程，距离中心点越近的部分，其权重越大，反之越小，而权值的分布，则符合正态分布，对于一个二维图像来说，则符合二维高斯分布。

高斯函数

一维高斯函数：

\[G(x) = \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{x^2}{2\sigma^2}} \]

二维高斯函数：

\[G(x,y) = \frac{1}{2\pi\sigma^2}e^{-\frac{x^2+y^2}{2\sigma^2}} \]

二维高斯函数的图像为：

技术图片

图像卷积

原理

为了下一步讨论高斯滤波的实现过程，我们首先要介绍一下图像处理中的卷积是如果工作和实现的。与信号处理中的卷积实现方法不同，在图像处理中的卷积主要通过模板矩阵与图像矩阵相乘，得到新的结果矩阵。下图详细演示了图像卷积的过程：

技术图片

在上图中，黄色标识的像素标识我们要卷积的目标像素，绿色为经过卷积核得出的结果像素值，其运算过程为：

\[0\times (-1) + 0\times 0 + 0\times 1+0\times (-2)+1\times 0+ 1\times 2 + 1\times(-1) + 1\times 0 + 1\times 1 = 3 \]

边缘补偿

在实际操作中，我们需要对图像的边缘进行补偿，以便卷积核可以对图像边缘部分的像素进行计算。我们常用的边缘补偿方法有以下几种：

补0：即补偿的像素值为0
复制：即复制边缘一圈的像素值

在本代码中，我们用补0的方式来进行处理。

代码

根据上文对图像处理卷积的理解，我们可以认为高斯卷积也是一个符合高斯分布的卷积核，在本实例中，我们使用的卷积核大小为\(3\times 3\)。

边缘补偿

如上文所述，在卷积处理开始之前，要对边缘进行补偿，我们这里使用全0补偿法，代码如下：

def padding(image,ksize):# 输入目标图像，卷积核大小

    h = image.shape[0] # 获取图像尺寸
    w = image.shape[1]
    c = image.shape[2]

    pad = ksize // 2 # 需要补偿的边缘大小

    out_p = np.zeros((h+2*pad,w+2*pad,c)) # 创造一个补偿后大小的全0矩阵

    out_copy = image.copy()

    out_p[pad:pad+h,pad:pad+w,0:c] = out_copy.astype(np.uint8) # 将原图像复制入目标图像中

    return out_p

高斯矩阵

下来是构建高斯卷积核，根据上图中的二维高斯函数公式，我们的代码为：

    # 高斯卷积核

    kernel = np.zeros((ksize,ksize)) # 创建一个卷积核大小的全0二维矩阵
    for x in range(-pad,-pad+ksize):
        for y in range(-pad,-pad+ksize):
            kernel[y+pad,x+pad] = np.exp(-(x**2+y**2)/(2*(sigma**2))) # 给卷积核内赋值
    kernel /= (sigma*np.sqrt(2*np.pi)) # 计算平均值
    kernel /=  kernel.sum() # 加总

要注意的是，高斯方程中有一个很重要的参数就是\(sigma\)值，在本示例中，如果未指定确定值，则\(sigma\)由下面的公式自动确定：

\[sigma = 0.3 \times ((ksize-1)\times 0.5-1) + 0.8 \]

其中ksize为高斯卷积核大小（本例中为\(3\)，即高斯卷积核为\(3\times 3\)）。

完整代码

最后，再根据图像中卷积的定义，对目标进行高斯卷积。

整个过程的完整代码如下：

import numpy as np 
from PIL import Image
from cv2 import cv2
import matplotlib.pyplot as plt
import math




def padding(image,ksize):

    h = image.shape[0]
    w = image.shape[1]
    c = image.shape[2]

    pad = ksize // 2 

    out_p = np.zeros((h+2*pad,w+2*pad,c))

    out_copy = image.copy()

    out_p[pad:pad+h,pad:pad+w,0:c] = out_copy.astype(np.uint8)

    return out_p





def gaussian(image,ksize,sigma):

    """
    1. padding
    2. 定义高斯滤波公式与卷积核
    3. 卷积过程

    高斯卷积卷积核是按照二维高斯分布规律产生的，公式为：

    G(x,y) = (1/(2*pi*(sigma)^2))*e^(-(x^2+y^2)/2*sigma^2)

    唯一的未知量是sigma，在未指定sigma的前提下，可以通过下列参考公式让程序自动选择合适的
    sigma值：

    sigma =  0.3 *((ksize-1)*0.5-1) + 0.8

    @ 如果mode为default，则返回abs值，否则返回unit8值
    """
        

    pad = ksize//2

    out_p = padding(image,ksize) # padding之后的图像
    # print(out_p)

    h = image.shape[0]
    w = image.shape[1]
    c = image.shape[2]

    # 高斯卷积核

    kernel = np.zeros((ksize,ksize))
    for x in range(-pad,-pad+ksize):
        for y in range(-pad,-pad+ksize):
            kernel[y+pad,x+pad] = np.exp(-(x**2+y**2)/(2*(sigma**2)))
    kernel /= (sigma*np.sqrt(2*np.pi))
    kernel /=  kernel.sum()

    # print(kernel)

    tmp = out_p.copy()

    # print(tmp)

    for y in range(h):
        for x in range(w):
            for z in range(c):

                out_p[pad+y,pad+x,z] = np.sum(kernel*tmp[y:y+ksize,x:x+ksize,z])


    out = out_p[pad:pad+h,pad:pad+w].astype(np.uint8)
    # print(out)

    return out
	
if __name__ == "__main__":

    path = "../bilder/lena.png"

    img = cv2.imread(path)
    gaussian_img = gaussian(img,3,0.8)
	cv2.imshow(‘Original Image‘,img)
	cv2.imshow(‘Gaussian Image‘,gaussian_img)
	cv2.waitkey()

结果

我们用图像处理中著名的"lena"图来进行结果展示。

技术图片

通过结果我们可以看出，高斯图像相比原图更加平滑，比如头发上的细节，经过高斯处理后，平滑了许多，如果想要更多的效果，可以通过增大sigma值来实现。

第二部分 Sobel算子

原理

Sobel算子是非常常用的边缘检测算子，那么，既然是边缘检测，首先就需要了解什么是边缘。

何为边缘？

我们人眼可以非常容易的分辨哪里是边缘哪里不是，但对于计算机来说，则没有那么容易。因为对于程序来说，一个灰度图像的本质是一个二维矩阵，矩阵中的每一个值对应的是图像每一个像素点的灰度值。那么，对于边缘来说，就是灰度变化最激烈的地方，换句话说，就是灰度梯度值最大的点，因为梯度就代表了变化程度。

Sobel算子也是一种加权算子，主要思想在于，各个像素对中心像素点的影响不是等价的，会随着距离的不同而不同。Sobel分别在x和y两个方向上对像素进行求导。

\[G(x) = \begin{bmatrix} -1&0&1\-2&0&2\-1&0&1 \end{bmatrix} * I \]

\[G(y) = \begin{bmatrix} -1&-2&-1\0&0&0\1&2&1 \end{bmatrix} * I \]

代码

在本例实现的代码中，我们定义了一个参数operator_type来确定是进行x还是y方向上的求导。如果operator_type = ‘horizontal‘则是x方向，operator_type = ‘Vertical‘则是y方向。不同求导方向对于不同方向上边缘的敏感性不同，具体的区别会在后面的结果清楚的呈现。

def sobel(image,operator_type,mode):


    if operator_type == "horizontal": # 定义求导方向

        sobel = np.array([[-1,-2,-1],[0,0,0],[1,2,1]])
    elif operator_type == "vertical":
        sobel = np.array([[-1,0,1],[-2,0,2],[-1,0,1]])
    else:
        print("type error") # 如果未定义，则输出错误提醒

    h = image.shape[0]
    w = image.shape[1]
    c = image.shape[2]
    out_p = padding(image,ksize=3)# padding，边缘补偿

    tmp = out_p.copy()

    for y in range(h):
        for x in range(w):
            for z in range(c):

                if mode == ‘default‘: # 为后面的canny算子进行准备，默认对结果像素取绝对值。

                    out_p[1+y,1+x,z] = abs(np.sum(sobel*tmp[y:y+3,x:x+3,z]))
                else:
                    out_p[1+y,1+x,z] = np.sum(sobel*tmp[y:y+3,x:x+3,z])

    out = out_p[1:1+h,1:1+w].astype(np.uint8)

    out = out.clip(0,255) # 裁剪像素，将结果像素值限定在0-255范围内。

    if mode == "default": # 设定阈值，当前阈值设定为70
        threshold = 70
        for i in range(h):
            for j in range(w):
                for z in range(c):
                    if out[i,j,z] < threshold:
                        out[i,j,z] = 0
    # print(out)

    return out