Python与FPGA——高斯滤波

前言

  从本节开始，高斯滤波和双边滤波都只进行Python的讲解，不涉及FPGA。还是一样，感兴趣的小伙伴可以自行实现。

一、高斯是谁？

  中心极限定理告诉我们，如果样本量最够大，样本的均值分布慢慢变成正太分布，那么图像的噪声也应该符合正太分布（高斯白噪声），所以就有了高斯滤波算法。至于高斯是谁并不重要，重要的是他提出了这一理论。公式如下。
$$\{\begin{array}{c}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}}(二维高斯分布)\end{array}$$

二、Python画出高斯分布

  一维高斯如下。

import numpy as np
import matplotlib.pyplot as plt
x = np.arange(-10, 10, 0.1)
sigma1 = 1
sigma2 = 2
sigma3 = 3
y1 = np.exp(-x ** 2 / (2 * sigma1 ** 2)) / ((2 * np.pi) ** 0.5 * sigma1)
y2 = np.exp(-x ** 2 / (2 * sigma2 ** 2)) / ((2 * np.pi) ** 0.5 * sigma2)
y3 = np.exp(-x ** 2 / (2 * sigma3 ** 2)) / ((2 * np.pi) ** 0.5 * sigma3)
plt.title("gauss curve")
plt.plot(x, y1)
plt.plot(x, y2)
plt.plot(x, y3)
plt.legend(["sigma1=1", "sigma2=2", "sigma3=3"])
plt.show()

  二维高斯如下。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D    #这里只是用Axes3D函数，所以只导入了Axes3D
fig = plt.figure()
ax = Axes3D(fig)
fig.add_axes(ax)
x = np.arange(-3, 3,0.05)     # 生成数据
y = np.arange(-3, 3, 0.05)
x,y=np.meshgrid(x,y)   #生成x,y轴数据
sigma = 1
z = np.exp(-(x ** 2 + y ** 2) / (2 * sigma ** 2)) / (2 * np.pi * sigma ** 2)
ax.plot_surface(x,y,z,rstride=1,cstride=1, cmap="rainbow")
ax.set_title("3D gauss")
ax.set_zlabel("z")
ax.set_xlabel("x")
ax.set_ylabel("y")
plt.show()

三、图像增噪

  为了方便，还是使用椒盐增噪，但是效果不是很好。

img = plt.imread("lenna.png")
gray = 0.299 * img[:, :, 0] + 0.587 * img[:, :, 1] + 0.114 * img[:, :, 2] 
gray = gray * 255#图像是[0-1]--->[0-255]
prob = 0.92#生成盐噪生的概率为1-0.92=0.08
noise = np.random.choice((0, 255), size=gray.shape, p=[prob, 1 - prob])#产生和图像维度一样的噪音矩阵
noise_img = np.where(noise == 0, gray, 255).astype(np.uint8)#生成噪音图像，椒（0黑色）盐噪声的盐（255白色）
fig = plt.figure(figsize=(10, 6))
ax = plt.subplot(1, 2, 1)
ax.set_title("salt image")
ax.set_xlabel("width")
ax.set_ylabel("height")
plt.imshow(gray, cmap="gray")
ax = plt.subplot(1, 2, 2)
ax.set_title("salt image")
ax.set_xlabel("width")
ax.set_ylabel("height")
plt.imshow(noise_img, cmap="gray")

四、Python高斯滤波

  $\sigma$不能负数，$\sigma$小于等于0的处理成0.3((ksize−1)/2−1)+0.8。

def gauss_kernel(kernel_size=3, sigma=1):
    k = kernel_size // 2
    if sigma == 0:
        sigma = ((kernel_size - 1) * 0.5 - 1) * 0.3 + 0.8
    x = np.linspace(-k, k, kernel_size)
    y = np.linspace(-k, k, kernel_size)
    x, y = np.meshgrid(x, y)
    z = np.exp(-(x ** 2 + y ** 2) / (2 * sigma ** 2)) / (2 * np.pi * sigma ** 2)
    return z

def guass_filter(image, kernel):
    h, w = image.shape
    filtered_image = np.zeros((h, w))
    n, n = kernel.shape
    m = int((n - 1) / 2)
    for i in range(m, h - m):
        for j in range(m, w - m):
           filtered_image[i, j] = np.sum(np.multiply(kernel, image[i - m: i + m + 1, j - m: j + m + 1]))
    return filtered_image.astype(np.uint8)
    
kernel = gauss_kernel(kernel_size=3)
guass_image = guass_filter(noise_img, kernel)
fig = plt.figure(figsize=(10, 6))
ax = plt.subplot(1, 2, 1)
ax.set_title("salt image")
ax.set_xlabel("width")
ax.set_ylabel("height")
plt.imshow(noise_img, cmap="gray")
ax = plt.subplot(1, 2, 2)
ax.set_title("gauss image")
ax.set_xlabel("width")
ax.set_ylabel("height")
plt.imshow(guass_image, cmap="gray")

总结

  从结果可以看到，高斯核对椒盐噪音无可奈何。如果你想让高斯滤波算法发挥功效，还是得使用有高斯噪音的图像。下一部分讲双边滤波，敬请期待。