图像平滑去噪之高斯滤波器

高斯滤波器是根据高斯函数来选择权值的线性平滑滤波器，对随机分布和服从正态分布的噪声有很好地滤除效果。本文从opencv内置的高斯滤波函数入手，深入介绍高斯滤波器的原理与实现。

一、高斯分布函数与高斯卷积核

高斯分布函数指的就是概率论中的正态分布的概率密度函数，均值μ=0时的一维形式和二维形式如下。其中σ为正态分布的标准偏差，其值决定了函数的衰减快慢。

从这两个公式不难看出，二维公式其实等于两个一维函数相乘。从概率论角度看，因为随机变量X，Y是相互独立的，那么他们的联合概率密度就等于边缘概率密度之积。这个特性是非常重要的，后面将再次提及。现在让我们先看一下高斯函数的图像分布与二维高斯卷积核的样子：

图像上，靠近原点的位置地势高，距离原点越远则地势越低。相应地，卷积核也是中心数值最大，并向四周减小，减小的幅度并不是随意的，而是要求整个卷积核近似高斯函数的图像。由于高斯滤波实质是一种加权平均滤波，为了实现平均，核还带有一个系数，例如上图中的十六分之一、八十四分之一。这些系数等于矩阵中所有数值之和的倒数。

二、高斯滤波

在实际滤波中，将当前像素作为核中心，利用卷积核对周围邻域像素作加权平均，其值作为当前像素的新值。这个过程实质上也是矩阵的卷积（矩阵卷积要求其中一个矩阵要旋转180度，在图像处理时，卷积核可以认为是已经过旋转的矩阵），因此才有卷积核这样的称呼。对于图像边界点，有多种处理方法，如对图像进行扩展、缩小高斯核等。OpenCV中提供了高斯滤波函数：

src为输入的图像，可以是单通道灰度图像，也可以是三通道RGB图像；ksize是卷积核大小，接受tuple类型；sigmaX是水平方向的标准差，sigmaY是垂直方向的标准差，这两个参数可以设置为0，函数内部会根据ksize自动计算标准差，标准差越大，核就越大，核变化趋势就越缓；borderType是边界的处理方式，采取默认值即可。以下是OpenCV-Python的滤波实例：

import cv2

img = cv2.imread('D:/programe/matlab/img/man.jpg')
gaussian = cv2.GaussianBlur(img,ksize=(5,5),sigmaX=0,sigmaY=0)
cv2.imshow('img',img)
cv2.imshow('gaussian',gaussian)
cv2.waitKey(0)

滤波前后对比如下图。因为高斯滤波是以空间距离来确定权重大小的，但并没有考虑用颜色距离来确定权重，这样就导致了高斯滤波在去除噪声的同时，也一定程度上模糊了边界。

三、高斯卷积核的构造方法

OpenCV本身就提供了高斯核的生成函数，不过它只能构建一维的核，但是根据上面我们提到过的一维乘积性质，我们可以通过两个一维核向量相乘得到一个二维核，该函数接口和代码如下：

#二维高斯卷积核
def GetGaussianKernel2D(ksize,sigma):
    kernelX = cv2.getGaussianKernel(ksize[0],sigma).reshape(ksize[0],1)#行向量
    kernelY = cv2.getGaussianKernel(ksize[1],sigma).reshape(1,ksize[1])#列向量
    kernel =  kernelX*kernelY
    return kernel
print(GetGaussianKernel2D((3,3),0.8))
#[[0.05711826 0.12475775 0.05711826]
# [0.12475775 0.27249597 0.12475775]
# [0.05711826 0.12475775 0.05711826]]

这里sigma是根据经验指定的，实际上更多地是让函数自动计算，标准差的计算公式为：

由于我们需要的是二维核，故可以设置两个标准差，分别表示X方向和Y方向，这两个标准差默认采用公式计算。最终获得的卷积核是小数形式，但是我们像素卷积需要整数，同时为了减小计算开销，我们希望核的值为较小的整数。因此，最常用的做法是将角点化为1，其它点相应乘以相同倍数再四舍五入或者直接去除小数部分，这样使得整个核的值较小，可以加快计算。由于标准差指定可能不同，因此还是会有很多相同大小但值不同的核存在，这并不需要惊讶，因为这些核都符合高斯分布，只是在实际滤波时效果有略微差异罢了。

#二维高斯卷积核
def GetGaussianKernel2D(ksize,sigmaX=None,sigmaY=None):
    if sigmaX is None:
        sigmaX = 0.3*((ksize[0]-1)*0.5-1)+0.8
    if sigmaY is None:
        sigmaY = 0.3 * ((ksize[1] - 1) * 0.5 - 1) + 0.8
    kernelX = cv2.getGaussianKernel(ksize[0],sigmaX).reshape(ksize[0],1)
    kernelY = cv2.getGaussianKernel(ksize[1],sigmaY).reshape(1,ksize[1])
    kernel =  kernelX*kernelY
    kernel =  kernelX*kernelY
    k = 1/kernel[0,0]
    kernel = np.int32(kernel*k+0.5)
    return kernel
print(GetGaussianKernel2D((3,3)))
#[[1 2 1]
# [2 5 2]
# [1 2 1]]

四、高斯卷积核的构造原理

我们也可以自己构造核，但首先我们要先知道高斯核与高斯函数之间的关系。以核中心为平面中心，则核的每个元素都对应一个坐标（x,y）如下图所示：

将每个元素对应的坐标代入到二维高斯函数中，得到的函数值作为该位置的值，再经归一化和整数化后就可以得到最终所要的卷积核。

首先构造高斯函数，Python中我没找到可以直接用的接口，所以就自定义算了，注意二维可以不用再写一遍公式，可以通过一维高斯获取：

#一维高斯函数
def Gaussian(x,sigma=1.0):
    result = 1/(np.sqrt(2*np.pi)*sigma)*np.exp(-(x/sigma)**2/2)
    return result

#二维高斯函数
def Gaussian2D(x,y,sigmaX=1.0,sigmaY=1.0):
    result = Gaussian(x,sigmaX)*Gaussian(y,sigmaY)
    return result

根据高斯核与高斯函数的关系，定义一维高斯核构造函数：

#一维高斯卷积核
def GaussianKernel(length,sigma=None):
    radius = length//2
    kernel = np.zeros(length,dtype=np.float32)
    if sigma is None:
        sigma = 0.3*((length-1)*0.5-1)+0.8
    for i in range(length):
        kernel.itemset(i,Gaussian(i-radius,sigma=sigma))
    #归一化
    sum = np.sum(kernel)
    kernel = kernel/sum
    return kernel

上面的函数就类似于OpenCV的getGaussianKernel。因此二维核也采用相同方法创建即可：

#二维高斯卷积核
def GaussianKernel2D(ksize,sigmaX=None,sigmaY=None):
    kernelX = GaussianKernel(ksize[0],sigmaX).reshape(ksize[0],1)
    kernelY = GaussianKernel(ksize[1],sigmaY).reshape(1,ksize[1])
    kernel =  kernelX*kernelY
    #整数化
    k = 1 / kernel[0, 0]
    kernel = np.int32(kernel * k + 0.5)
    return kernel

五、卷积可分离性

在高斯卷积应用中，可能会用到较大的卷积核，此时程序性能就不太令人满意了。卷积可分离性就是解决这种问题的一种优化手段。对于高斯卷积来说，我们可以将图像先与一维水平卷积核相卷积，再与一维垂直方向的卷积核相卷积，所得结果跟直接用二维核来卷积相同，这就是高斯卷积核的可分离性。利用该分离性，卷积复杂度大大降低，并且我们也不需要生成二维卷积核，而只要生成一维核即可，这也减少了计算量。在实际应用中，甚至还可以直接存储常用维数的高斯核数值来提升性能。

下面谈一谈卷积可分离性的原理，不感兴趣的读者可直接阅读下一节。

假设A为n维列向量，B为n维行向量，则有A*B（卷积）=B*A=A·B。例如n=3时,如下图所示：

A对B的卷积（A作为核）需要注意先将A向量旋转180度，再做加权求和操作，过程如下：

前面说过高斯核Kernel2D=KernelX·KernelY，那么卷积过程就可以表示为：Dst=Src*Kernel2D=(Src*KernelX)*KernelY=(Src*KernelY)*KernelX。一般来说，无论对于何种卷积滤波，只要其卷积核可以拆解为两个行向量与列向量的乘积，那么就算卷积可分离。

六、优化的高斯滤波实现

这一节综合了上面所有代码，并且添加了优化的自定义高斯滤波函数：

import numpy as np
import cv2


#一维高斯函数
def Gaussian(x,sigma=1.0):
    result = 1/(np.sqrt(2*np.pi)*sigma)*np.exp(-(x/sigma)**2/2)
    return result

#二维高斯函数
def Gaussian2D(x,y,sigmaX=1.0,sigmaY=1.0):
    result = Gaussian(x,sigmaX)*Gaussian(y,sigmaY)
    return result

#一维高斯卷积核
def GaussianKernel(length,sigma=None):
    radius = length//2
    kernel = np.zeros(length,dtype=np.float32)
    if sigma is None:
        sigma = 0.3*((length-1)*0.5-1)+0.8
    for i in range(length):
        kernel.itemset(i,Gaussian(i-radius,sigma=sigma))
    #归一化
    sum = np.sum(kernel)
    kernel = kernel/sum
    return kernel

#二维高斯卷积核
def GaussianKernel2D(ksize,sigmaX=None,sigmaY=None):
    kernelX = GaussianKernel(ksize[0],sigmaX).reshape(ksize[0],1)
    kernelY = GaussianKernel(ksize[1],sigmaY).reshape(1,ksize[1])
    kernel =  kernelX*kernelY
    k = 1 / kernel[0, 0]
    kernel = np.int32(kernel * k + 0.5)
    return kernel

#二维高斯卷积核
def GetGaussianKernel2D(ksize,sigmaX=None,sigmaY=None):
    if sigmaX is None:
        sigmaX = 0.3*((ksize[0]-1)*0.5-1)+0.8
    if sigmaY is None:
        sigmaY = 0.3 * ((ksize[1] - 1) * 0.5 - 1) + 0.8
    kernelX = cv2.getGaussianKernel(ksize[0],sigmaX).reshape(ksize[0],1)
    kernelY = cv2.getGaussianKernel(ksize[1],sigmaY).reshape(1,ksize[1])
    kernel =  kernelX*kernelY
    k = 1/kernel[0,0]
    kernel = np.int32(kernel*k+0.5)
    return kernel

#优化的高斯滤波
def GaussianFilter(src,ksize):
    #卷积分离
    kernelX = GaussianKernel(ksize[0])
    kernelX = np.resize(kernelX,(ksize[0],1))
    kernelY = GaussianKernel(ksize[1])
    kernelX = np.resize(kernelX, (1,ksize[1]))
    dst = cv2.filter2D(src,-1,kernelX)
    dst = cv2.filter2D(dst,-1,kernelY)
    return dst


if __name__ == '__main__':
    src=cv2.imread('man.jpg')
    dst = GaussianFilter(src,(5,5))
    cv2.imshow('src',src)
    cv2.imshow('dst',dst)
    cv2.waitKey(0)

posted @ 2019-10-06 18:40 KenSporger 阅读(21028) 评论(0) 编辑收藏举报

刷新页面返回顶部

KenSporger