第3章 Python 数字图像处理(DIP) - 灰度变换与空间滤波14 - 平滑低通滤波器 -高斯滤波器核的生成方法

最新推荐文章于 2024-07-07 22:02:29 发布

jasneik

最新推荐文章于 2024-07-07 22:02:29 发布

阅读量896

点赞数 1

分类专栏：数字图像处理(DIP)--Python # 第3章 - 灰度变换与空间滤波文章标签：边缘检测图像处理 numpy opencv 高斯拉普拉斯算子

本文链接：https://blog.csdn.net/jasneik/article/details/113767260

版权

数字图像处理(DIP)--Python 同时被 2 个专栏收录

57 篇文章 128 订阅

订阅专栏

第3章 - 灰度变换与空间滤波

18 篇文章 6 订阅

订阅专栏

平滑（低通）空间滤波器

平滑（也称平均）空间滤波器用于降低灰度的急剧过渡

在图像重取样之前平滑图像以减少混淆
用于减少图像中无关细节
平滑因灰度级数量不足导致的图像中的伪轮廓
平滑核与一幅图像的卷积会模糊图像

低通高斯滤波器核

高斯核是唯一可分离的圆对称（也称各向同性，意味它们的响应与方向无关）核。
$e^{-\frac{s^2 + t^2}{2\sigma^2}} \tag{3.45}$
令 $r = [s^2 + t^2]^{1/2}$ ，上式可以重写为
$e^{-\frac{r^2}{2\sigma^2}} \tag{3.46}$

高斯分布函数
$\frac{1}{\sqrt{2\pi\sigma}}e^{-\frac{(x-x_0)^2}{2\sigma^2}}$
$\frac{1}{2\pi\sigma^2}e^{-\frac{(x-x_0)^2+(y-y_0)^2}{2\sigma^2}}$

3X3高斯核：

$\frac{1}{4.897} *\begin{bmatrix} 0.3679 & 0.6065 & 0.3679 \\ 0.6065 & 1 & 0.6065\\ 0.3679 & 0.6065 & 0.3679 \end{bmatrix}$

下面的近似高斯核（效果很接近）

$\frac{1}{16} *\begin{bmatrix} 1& 2& 1 \\ 2& 4& 2 \\ 1& 2& 1\end{bmatrix}$

$\frac{1}{273} * \begin{bmatrix} 1& 4& 7& 4& 1 \\ 4& 16& 26& 16& 4 \\ 7& 26& 41& 26& 7 \\ 4& 16& 26& 16& 4 \\ 1& 4& 7&4&1\end{bmatrix}$

下面是两种方式构建低通高斯滤波器核

def gauss_seperate_kernel(kernel_size, sigma=1):
    """
    create gaussian kernel use separabal vector, because gaussian kernel is symmetric and separable.
    param: kernel_size: input, the kernel size you want to create, normally is uneven number, such as 1, 3, 5, 7
    param: sigma: input, sigma of the gaussian fuction, default is 1
    return a desired size 2D gaussian kernel
    """
    if sigma == 0:
        sigma = ((kernel_size - 1) * 0.5 - 1) * 0.3 + 0.8
    m = kernel_size[0]
    n = kernel_size[1]
    x_m = np.arange(m) - m // 2
    y_m = np.exp(-(x_m **2 )/ (2 * sigma**2)).reshape(-1, 1)
    
    x_n = np.arange(n) - n // 2
    y_n = np.exp(-(x_n ** 2)/ (2 * sigma**2)).reshape(-1, 1)
    kernel = y_m * y_n.T
    return kernel

# gaussian kernel
gauss_seperate_kernel((5, 5)).round(3)

array([[0.018, 0.082, 0.135, 0.082, 0.018],
       [0.082, 0.368, 0.607, 0.368, 0.082],
       [0.135, 0.607, 1.   , 0.607, 0.135],
       [0.082, 0.368, 0.607, 0.368, 0.082],
       [0.018, 0.082, 0.135, 0.082, 0.018]])

def gauss_kernel(kernel_size, sigma=1):
    """
    create gaussian kernel use meshgrid, gaussian is circle symmetric.
    param: kernel_size: input, the kernel size you want to create, normally is uneven number, such as 1, 3, 5, 7
    param: sigma: input, sigma of the gaussian fuction, default is 1
    return a desired size 2D gaussian kernel
    """
#     X = np.arange(-kernel_size // 2 +1, kernel_size//2 +1)
#     Y = np.arange(-kernel_size // 2 +1, kernel_size//2 +1)
    
    m = kernel_size[0]
    n = kernel_size[1]
    X = np.arange(- m//2 + 1, m//2+1)
    Y = np.arange(- n//2 + 1, n//2+1)
    x, y = np.meshgrid(X, Y)
    kernel = np.exp(- ((x)**2 + (y)**2)/ (2 * sigma**2))
    return kernel

# gaussian kernel
gauss_kernel((5, 5)).round(3)

array([[0.018, 0.082, 0.135, 0.082, 0.018],
       [0.082, 0.368, 0.607, 0.368, 0.082],
       [0.135, 0.607, 1.   , 0.607, 0.135],
       [0.082, 0.368, 0.607, 0.368, 0.082],
       [0.018, 0.082, 0.135, 0.082, 0.018]])

def visualize_show_annot(img_show, img_annot, ax):
    """
    add annotation to the image, values of each pixel
    param: img: input image
    param: ax: axes of the matplotlib
    """
    height, width = img_annot.shape
    img_show = img_show[:height, :width]
    ax.imshow(img_show, cmap='gray', vmin=0, vmax=255)
    thresh = 10 #img_show.max()/2.5
    for x in range(height):
        for y in range(width):
            ax.annotate(str(round(img_annot[x][y],2)), xy=(y,x),
                        horizontalalignment='center',
                        verticalalignment='center',
                        color='white' if img_annot[x][y]>thresh else 'black')

# 显示距离
m = 9
S = np.arange(- (m - 1) // 2, (m - 1) // 2 + 1)
T = S
s, t = np.meshgrid(S, T)
r = np.power((s**2 + t**2), 0.5)

height, width = m, m
img_show = np.ones([height, width], dtype=np.uint8) * 250
fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot(1, 1, 1)
visualize_show_annot(img_show, r, ax)
plt.xticks([0, m-1], labels=["\sqrt{2}(m - 1)/2", "\sqrt{2}(m - 1)/2"]) # 标签显示不了\sqrt{2}，所以简写(m-1)/2
ax.set_yticks([0, m-1])
ax.set_yticklabels(["\sqrt{2}(m - 1)/2", "\sqrt{2}(m - 1)/2"])
plt.show()

在这里插入图片描述

# 显示高斯核滤波器的圆
import numpy as np
import mpl_toolkits.mplot3d
from matplotlib import pyplot as plt
from matplotlib import cm

x, y = np.mgrid[-3:3:200j,-3:3:200j]
z = np.exp(- ((x)**2 + (y)**2)/ (2))

fig = plt.figure(figsize=(8,8))
ax = fig.gca()
ax.contour(x, y, z, cmap=cm.ocean)
plt.show()

在这里插入图片描述

# 二维高斯核函数
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from matplotlib import pyplot as plt

k = 3
sigma = 1
X = np.linspace(-k, k, 100)
Y = np.linspace(-k, k, 100)
x, y = np.meshgrid(X, Y)
z = np.exp(- ((x)**2 + (y)**2)/ (2 * sigma**2))

fig = plt.figure(figsize=(8, 8))
ax = Axes3D(fig)
ax.set_xticks([])
ax.set_yticks([])
ax.set_zticks([])
ax.set_xlabel('s', fontsize=14)
ax.set_ylabel('t', fontsize=14)
ax.set_zlabel('G(s,t)', fontsize=14)
ax.plot_surface(x, y, z, antialiased=True, shade=False, alpha=0.5, cmap=cm.terrain)
ax.contour(x, y, z, [0.3, 0.6], colors='k')
ax.view_init(30, 60)
plt.show()

在这里插入图片描述

两个一维高斯函数 $f$ 和 $g$ 和乘积( $\times$ )和卷积( $\star$ )的均值与标准差

	$f$	$g$	$f\times{g}$	$f\star{g}$
均值	$m_f$	$m_g$	$m_{f\times{g}}=\frac{m_{f} \sigma_{g}^2 + m_{g} \sigma_{f}^2}{\sigma_{f}^2 + \sigma_{g}^2}$	$m_{f\star{g}} = m_{f} + m_{g}$
标准差	$\sigma_{f}$	$\sigma_{g}$	$\sigma_{f\times{g}}=\sqrt{\frac{\sigma_{f}^2 \sigma_{g}^2}{\sigma_{f}^2 + \sigma_{g}^2}}$	$\sigma_{f\star{g}} = \sigma_{f}^2 + \sigma_{g}^2$

# 高斯核生成函数 之前没有怎么看书的时候写的，不太正确
def creat_gauss_kernel(kernel_size=3, sigma=1, k=1):
    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)
    x0 = 0
    y0 = 0
    gauss = 1/(2*np.pi*sigma**2) * np.exp(- ((x -x0)**2 + (y - y0)**2)/ (2 * sigma**2))
    
##绘图功能
#     fig = plt.figure(figsize=(8,6))
#     ax = fig.gca(projection='3d')
#     ax.plot_surface(x,y,gauss,cmap=plt.cm.ocean)
#     plt.show()
    return gauss

#不同的大小，不同sigma的高斯核
img_gray = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH03/Fig0333(a)(test_pattern_blurring_orig).tif', 0)
kernel_21 = gauss_kernel((21, 21), sigma=3.5)
kernel_21 = kernel_21 / kernel_21.sum()
kernel_43 = gauss_kernel((43, 43), sigma=7)
kernel_43 = kernel_43 / kernel_43.sum()
# img_21 = cv2.filter2D(img_gray, ddepth= -1, kernel=kernel_21)
# img_43 = cv2.filter2D(img_gray, ddepth= -1, kernel=kernel_43)
img_21 = my_conv2d(img_gray, kernel=kernel_21)
img_43 = my_conv2d(img_gray, kernel=kernel_43)

plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1), plt.imshow(img_gray, cmap='gray', vmin=0, vmax=255), plt.title(f'Original')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 2), plt.imshow(img_21, cmap='gray', vmin=0, vmax=255), plt.title(f'Gaussian Kernel 21x21')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 3), plt.imshow(img_43, cmap='gray', vmin=0, vmax=255), plt.title(f'Gaussian Kernel 43x43')
plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

#不同的大小，不同sigma的高斯核，用可分离核运算的话，高斯核不用归一化，但计算后归一化后再乘以255
img_gray = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH03/Fig0333(a)(test_pattern_blurring_orig).tif', 0)
kernel_21 = gauss_kernel((21, 21), sigma=3.5)
# kernel_21 = kernel_21 / kernel_21.sum()
kernel_43 = gauss_kernel((43, 43), sigma=7)
# kernel_43 = kernel_43 / kernel_43.sum()
img_21 = separate_kernel_conv2D(img_gray, kernel=kernel_21)
img_21 = np.uint8(normalize(img_21) * 255)
img_43 = separate_kernel_conv2D(img_gray, kernel=kernel_43)
img_43 =  np.uint8(normalize(img_43) * 255)
plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1), plt.imshow(img_gray, cmap='gray', vmin=0, vmax=255), plt.title(f'Original')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 2), plt.imshow(img_21, cmap='gray', vmin=0, vmax=255), plt.title(f'Gaussian Kernel 21x21')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 3), plt.imshow(img_43, cmap='gray', vmin=0, vmax=255), plt.title(f'Gaussian Kernel 43x43')
plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

#大于[6sigma] x [6sigma]的高斯核并无益处，由于opencv返回的图像都是uint8的，从差来看，还是有点区别，但是两者的结果视觉上看差不多
#用自己实现的卷积，返回的是float，结果的最大值为0.948，比opencv的效果好，就是执行的效率不高，opencv计算是相关性，而不是真的卷积运算
img_gray = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH03/Fig0333(a)(test_pattern_blurring_orig).tif', 0)

kernel_43 = gauss_kernel((43, 43), sigma=7)
kernel_43 = kernel_43 / kernel_43.sum()
# img_43 = cv2.filter2D(img_gray, ddepth= -1, kernel=kernel_43)
img_43 = my_conv2d(img_gray, kernel=kernel_43)

kernel_85 = gauss_kernel((85, 85), sigma=7)
kernel_85 = kernel_85 / kernel_85.sum()
# img_85 = cv2.filter2D(img_gray, ddepth= -1, kernel=kernel_85)
img_85 = my_conv2d(img_gray, kernel=kernel_85)

img_diff = img_85 - img_43
plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1), plt.imshow(img_43, cmap='gray', vmin=0, vmax=255), plt.title(f'Gaussian Kernel 43x43')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 2), plt.imshow(img_85, cmap='gray', vmin=0, vmax=255), plt.title(f'Gaussian Kernel 85x85')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 3), plt.imshow(img_diff, cmap='gray', vmin=0, vmax=255), plt.title(f'Difference')
plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

def down_sample(image):
    height, width = image.shape[:2]
    dst = np.zeros([height//2, width//2])
    dst = image[::2, ::2]
    return dst

#高斯核与盒式平滑特性的比较，图像原尺寸为1024，但报内错误，做降采样，可以运行了。
#矩形区域为768x128，线的宽度为10，
height, width = 1024, 1024
img_ori = np.zeros([height, width], np.uint8)
shift = 120
img_ori[shift:shift+768, shift:shift+128] = 255
img_ori[610:620, 110:258] = 255   # 矩形背后的横线
img_ori[610:620, 300:500] = 255   #横线
img_ori[200:620, 500:510] = 255   #竖线
img_ori[200:210, 510:660] = 255   #横线
img_ori[210:620, 650:660] = 255   #竖线
img_ori[610:620, 660:]    = 255   #横线

# 做降采样缩小图像可以不报错了
img_ori = down_sample(img_ori)
# 盒式滤波器核 35x35，书上是71x71
box_71 = np.ones([71, 71], np.float)
box_71 = box_71 / box_71.sum()
# 高斯滤波器核71x71，书上要求是151x151，但会报内存错误
gauss_151 = gauss_kernel((71, 71), sigma=25)
gauss_151 = gauss_151 / gauss_151.sum()
img_box_71 = my_conv2d(img_ori, kernel=box_71)
img_gauss_151 = my_conv2d(img_ori, kernel=gauss_151)

# # # opencv处理线都不见了，不是真的做卷积运算
# # img_box_71 = cv2.filter2D(img_ori, ddepth=-1, kernel=box_71)
# # img_gauss_151 = cv2.filter2D(img_ori, ddepth=-1, kernel=gauss_151)

plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1), plt.imshow(img_ori, cmap='gray', vmin=0, vmax=255), plt.title(f'Original')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 2), plt.imshow(img_box_71, cmap='gray', vmin=0, vmax=255), plt.title(f'Box Kernel 71x71')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 3), plt.imshow(img_gauss_151, cmap='gray', vmin=0, vmax=255), plt.title(f'gaussian Kernel 151x151')
plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

def nomalize(img):
    return (img - img.min()) / (img.max() - img.min())

#高斯核与盒式平滑特性的比较，图像原尺寸为1024，用可分离核计算1024，也不会报错的。
#矩形区域为768x128，线的宽度为10，
height, width = 1024, 1024
img_ori = np.zeros([height, width], np.uint8)
shift = 120
img_ori[shift:shift+768, shift:shift+128] = 255
img_ori[610:620, 110:258] = 255   # 矩形背后的横线
img_ori[610:620, 300:500] = 255   #横线
img_ori[200:620, 500:510] = 255   #竖线
img_ori[200:210, 510:660] = 255   #横线
img_ori[210:620, 650:660] = 255   #竖线
img_ori[610:620, 660:]    = 255   #横线

# 做降采样缩小图像可以不报错了
# img_ori = down_sample(img_ori)
# 盒式滤波器核 35x35，书上是71x71
box_71 = np.ones([71, 71], np.float)
box_71 = box_71 / box_71.sum()
# 高斯滤波器核71x71，书上要求是151x151，但会报内存错误
gauss_151 = gauss_kernel((71, 71), sigma=25)
gauss_151 = gauss_151 / gauss_151.sum()
img_box_71 = separate_kernel_conv2D(img_ori, kernel=box_71)
img_box_71 = np.uint8(nomalize(img_box_71) * 255)
img_gauss_151 = separate_kernel_conv2D(img_ori, kernel=gauss_151)
img_gauss_151 = np.uint8(nomalize(img_gauss_151) * 255)
# # # opencv处理线都不见了，不是真的做卷积运算
# # img_box_71 = cv2.filter2D(img_ori, ddepth=-1, kernel=box_71)
# # img_gauss_151 = cv2.filter2D(img_ori, ddepth=-1, kernel=gauss_151)

plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1), plt.imshow(img_ori, cmap='gray', vmin=0, vmax=255), plt.title(f'Original')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 2), plt.imshow(img_box_71, cmap='gray', vmin=0, vmax=255), plt.title(f'Box Kernel 71x71')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 3), plt.imshow(img_gauss_151, cmap='gray', vmin=0, vmax=255), plt.title(f'gaussian Kernel 151x151')
plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

#相同的核，不同的边填充对结果的边的影响
img_gray = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH03/Fig0333(a)(test_pattern_blurring_orig).tif', 0)

# kernel size is 111x111, sigma=31，核太大了，计算会消耗很大的内存
kernel_111 = guass_kernel((111, 111), sigma=15)
kernel_111 = kernel_111 / kernel_111.sum()

img_constant = my_conv2d(img_gray, kernel=kernel_111, mode='constant')
img_reflect = my_conv2d(img_gray, kernel=kernel_111, mode='reflect')
img_copy = my_conv2d(img_gray, kernel=kernel_111, mode='edge')

plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1), plt.imshow(img_constant, cmap='gray', vmin=0, vmax=255), plt.title(f'0 padding')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 2), plt.imshow(img_reflect, cmap='gray', vmin=0, vmax=255), plt.title(f'Reflect padding')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 3), plt.imshow(img_copy, cmap='gray', vmin=0, vmax=255), plt.title(f'Copy padding')
plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

def up_sample_2d(img):
    """
    up sample 2d image, if your image is RGB, you could up sample each channel, then use np.dstack to merge a RGB image
    param: input img: it's a 2d gray image
    return a 2x up sample image
    """
    height, width = img.shape[:2]
    temp = np.zeros([height*2, width*2])
    temp[::2, ::2] = img
    temp[1::2, 1::2] = img
    temp[0::2, 1::2] = img
    temp[1::2, 0::2] = img
    return temp

#核和图像大小对平滑性能的影响
img_gray = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH03/Fig0333(a)(test_pattern_blurring_orig).tif', 0)

# 原图像是500x500，结果两次上采样，变是2000x2000
img_up = up_sample_2d(img_gray)
# img_up = up_sample_2d(img_up)

# 本想要升到2000x2000再做卷积，但报内存错误，只好试1000x1000的了
kernel_43 = guass_kernel((43, 43), sigma=7)
kernel_43 = kernel_43 / kernel_43.sum()
img_43 = my_conv2d(img_up, kernel=kernel_43, mode='reflect')

# kernel size is 111x111, sigma=31，核太大了，计算会消耗很大的内存
kernel_111 = guass_kernel((111, 111), sigma=15)
kernel_111 = kernel_111 / kernel_111.sum()
img_111 = my_conv2d(img_gray, kernel=kernel_111, mode='reflect')

plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1), plt.imshow(img_gray, cmap='gray', vmin=0, vmax=255), plt.title(f'0 padding')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 2), plt.imshow(img_43, cmap='gray', vmin=0, vmax=255), plt.title(f'Reflect padding Kernel 43x43')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 3), plt.imshow(img_111, cmap='gray', vmin=0, vmax=255), plt.title(f'Reflect padding Kernel 111x111')
plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

#使用低通滤波和阈值处理提取区域，根据自己的要求要选择核的大小
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH03/Fig0334(a)(hubble-original).tif', 0)
img_ori = img_ori / 255.

kernel_43 = gauss_kernel((43, 43), sigma=7)
kernel_43 = kernel_43 / kernel_43.sum()
img_43 = my_conv2d(img_ori, kernel=kernel_43, mode='reflect')

thred = 0.42
img_thred = np.where(img_43 <= thred, img_43, 1)
img_thred = np.where(img_43 > thred, img_thred, 0)

plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1), plt.imshow(img_ori, cmap='gray', vmin=0, vmax=1), plt.title(f'Original')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 2), plt.imshow(img_43, cmap='gray', vmin=0, vmax=1), plt.title(f'Gaussian filter')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 3), plt.imshow(img_thred, cmap='gray', vmin=0, vmax=1), plt.title(f'Binary')
plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

阴影校正（也称平坦场校正）

def filter_2d(image, kernel, mode='constant'):
    """
    convolution 2d image with kernel
    :param image: input image
    :param kernel: input kernel
    :return: image after convolution
    """
    
    img_h = image.shape[0]
    img_w = image.shape[1]

    m = kernel.shape[0]
    n = kernel.shape[1]

    # padding
    padding_h = int((m -1)/2)
    padding_w = int((n -1)/2)
    
    image_pad = np.pad(image, [(padding_h, padding_h), (padding_w, padding_w)], mode=mode)
#     image_pad = np.zeros((image.shape[0]+padding_h*2, image.shape[1]+padding_w*2), np.uint8)
#     image_pad[padding_h:padding_h+img_h, padding_w:padding_w+img_w] = image
    
    image_convol = image.copy()
    for i in range(padding_h, img_h + padding_h):
        for j in range(padding_w, img_w + padding_w):
            temp = np.sum(image_pad[i-padding_h:i+padding_h+1, j-padding_w:j+padding_w+1] * kernel)
            image_convol[i - padding_h][j - padding_w] = temp # 1/(m * n) * temp

    return image_convol

def get_check(height, width, check_size=(5, 5), lower=130, upper=255):
    """
    create check pattern image
    height: input, height of the image you want
    width: input, width of the image you want
    check_size: the check size you want, default is 5x5
    lower: dark color of the check, default is 130, which is dark gray, 0 is black, 255 is white
    upper: light color of the check, default is 255, which is white, 0 is black
    return uint8[0, 255] grayscale check pattern image
    """
    m, n = check_size
    black = np.zeros((m, n), np.uint8)
    white = np.zeros((m, n), np.uint8)
    black[:] = lower  # dark
    white[:] = upper  # white

    black_white = np.concatenate([black, white], axis=1)
    white_black = np.concatenate([white, black], axis=1)
    black_white_black_white = np.vstack((black_white, white_black))

    tile_times_h = int(np.ceil(height / m / 2))
    tile_times_w = int(np.ceil(width / n / 2))

    img_temp = np.tile(black_white_black_white, (tile_times_h, tile_times_w))
    img_dst = np.zeros([height, width])
    img_dst = img_temp[:height, :width]
    
    return img_dst

def get_shape_pattern(img_ori):
    """
    create a dark mask shape of the image
    img_ori: input image, which only use the image shape here to create the mask
    return [0, 1] float image
    """
    height, width = img_ori.shape
    img_dst = np.zeros([height, width])
    for i in range(height):
        temp = np.linspace(0, i/width, height)
        img_dst[i, :] = temp
    img_dst = normalize(img_dst)
    return img_dst

#使用低通滤波校正阴影，大核太慢了，还是用opencv的快
#创建棋盘格
height, width = 2048, 2048
img_ori = get_check(height, width, check_size=(128, 128), lower=0)

#得到棋盘格图像对应的掩膜
img_shape = get_shape_pattern(img_ori)

#生成带遮照的图像
img_ori = img_ori * img_shape

#生成513x513的高斯滤波器核
kernel_513 = guass_kernel((511, 511), sigma=128)
kernel_513 = kernel_513 / kernel_513.sum()

#卷积，别用注释掉的，按书上说的得2个小时，我也没有试过
# img_513 = filter_2d(img_ori, kernel=kernel_513)
img_513 = cv2.filter2D(img_ori, ddepth=-1, kernel=kernel_513)

#得到校正后的图像
img_dst = img_ori / (img_513 + 1e-5)

plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1), plt.imshow(img_ori, cmap='gray', vmin=0, vmax=255), plt.title(f'Original')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 2), plt.imshow(img_513, cmap='gray', vmin=0, vmax=1), plt.title(f'After Guassian Filter')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 3), plt.imshow(img_dst, cmap='gray', vmin=0, vmax=1), plt.title(f'Calibrated')
plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

# 多种滤波器核
img = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH03/Fig0333(a)(test_pattern_blurring_orig).tif')
img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# 盒式滤波器核
kernel_normal = np.array([[1, 1, 1],
                         [1, 1, 1],
                         [1, 1, 1]], dtype='float32')
kernel_normal = kernel_normal / kernel_normal.sum()

# 模糊算子
kernal_blur = np.array((
        [0.0625, 0.125, 0.0625],
        [0.125, 0.25, 0.125],
        [0.0625, 0.125, 0.0625]), dtype="float32")
kernal_blur = kernal_blur / kernal_blur.sum()

# sobel算子，sobel内核用于仅显示特定方向上相邻像素值的差异，分为left sobel、right sobel（检测梯度的水平变化）、top sobel、buttom sobel（检测梯度的垂直变化）。
kernel_sobel = np.array((
                    [-1,-2,-1],
                    [0,0,0],
                    [1,2,1]), np.int8)
# 浮雕算子，通过强调像素的差在给定方向的Givens深度的错觉。在这种情况下，沿着从左上到右下的直线的方向。
kernel_emboss = np.array((
                    [-2,-1,0],
                    [-1,1,1],
                    [0,-1,2]), np.int8)

# 轮廓内核（也称为“边缘”的内核）用于突出显示的像素值大的差异。具有接近相同强度的相邻像素旁边的像素在新图像中将显示为黑色，而与强烈不同的相邻像素相邻的像素将显示为白色。
kernel_outline = np.array((
                    [-1,-1,-1],
                    [-1,8,-1],
                    [-1,-1,-1]), np.int8)

# 锐化算子，锐化内核强调在相邻的像素值的差异，这使图像看起来更生动
kernel_shapen = np.array((
                    [0,-1,0],
                    [-1,5,-1],
                    [0,-1,0]), np.int8)

# 拉普拉期算子，用于边缘检对于检测图像中的模糊也非常有用
kernel_laplacian = np.array((
                    [0,1,0],
                    [1,-4,1],
                    [0,1,0]), np.int8)

# 分身算子，就是原图
kernel_identity = np.array((
                    [0,0,0],
                    [0,1,0],
                    [0,0,0]), np.int8)

# 网上找的高斯核
kernel_gaussianfake = np.array([[1, 2, 1],
                                [2, 4, 2],
                                [1, 2, 1]])/16
# DIP书上给了的高斯核
kernel_gaussian = np.array([[0.3679, 0.6065, 0.3679],
                            [0.6065, 1, 0.6065],
                            [0.3679, 0.6065, 0.3679]], np.float32) / 4.8976

# 自己通过高期函数算出的高斯核
kernel_creategaussian = guass_kernel((3, 3))
kernel_creategaussian = kernel_creategaussian / kernel_creategaussian.sum()

img_nomal = cv2.filter2D(img_gray, ddepth= -1, kernel=kernel_normal)
imgkernal_blur = cv2.filter2D(img_gray, -1, kernal_blur)
imgkernel_sobel = cv2.filter2D(img_gray, -1, kernel_sobel)
imgkernel_emboss = cv2.filter2D(img_gray, -1, kernel_emboss)
imgkernel_outline = cv2.filter2D(img_gray, -1, kernel_outline)
imgkernel_shapen = cv2.filter2D(img_gray, -1, kernel_shapen)
imgkernel_laplacian = cv2.filter2D(img_gray, -1, kernel_laplacian)
imgkernel_identity = cv2.filter2D(img_gray, -1, kernel_identity)
imgkernel_gaussianfake = cv2.filter2D(img_gray, -1, kernel_gaussianfake)
imgkernel_gaussian = cv2.filter2D(img_gray, -1, kernel_gaussian)
imgkernel_creategaussian = cv2.filter2D(img_gray, -1, kernel_creategaussian)
img_list = ['img_gray', 'img_nomal', 'imgkernal_blur', 'imgkernel_identity', 'imgkernel_gaussian', 'imgkernel_creategaussian', 
            'imgkernel_shapen', 'imgkernel_gaussianfake', 'imgkernel_sobel', 'imgkernel_emboss', 'imgkernel_outline', 
            'imgkernel_laplacian', ]

fig = plt.figure(figsize=(15, 20))

for i, img in enumerate(img_list):
    ax = fig.add_subplot(4, 3, i+1)
    ax.imshow(locals()[img], cmap='gray')
    ax.set_title(img.split('_')[-1])
    
plt.tight_layout()
plt.show()

在这里插入图片描述

# 拉普拉斯算子锐化
laplacian_img = np.uint8(normalize(img_gray - imgkernel_laplacian) * 255)

plt.figure(figsize=(18, 6))
plt.subplot(1, 3, 1)
plt.imshow(img_gray, cmap='gray')
plt.title('Original')

plt.subplot(1, 3, 2)
plt.imshow(laplacian_img, cmap='gray')
plt.title('Laplacian')

img_diff = img_gray - laplacian_img
plt.subplot(1, 3, 3)
plt.imshow(img_diff, cmap='gray')
plt.title('Diff')
plt.tight_layout()
plt.show()

在这里插入图片描述

统计排序（非线性）滤波器

统计排序滤波器是非线性空间滤波器，其响应基于滤波器所包含区域内的像素的排序。平滑是将中心像素的值替代为由排序结果确定的值来实现的。

中值滤波器
中值滤波器对某些类型的随机噪声提供了优秀的降噪能力
与类似大小的线性平滑滤波器相对，中值滤波器对图像的模糊程度要小得多。
对椒盐噪声尤其有效
邻域内的中值

def median_filter(img, kernel_size, mode="constant"):
    """
    median filter
    param: img: input image for denoising
    param: kernel_size: kernel size of the median filter block
    return: image after median filter
    """
    img_h = img.shape[0]
    img_w = img.shape[1]
    
    m = kernel_size[0]
    n = kernel_size[1]
    
    pad_h = int((m - 1)/2)
    pad_w = int((n - 1)/2)
    
    # 边界通过0灰度值填充扩展
    img_pad = np.pad(img, [(pad_h, pad_h), (pad_w, pad_w)], mode=mode)
    
    img_result = np.zeros(img.shape)
    for i in range(img_result.shape[0]):
        for j in range(img_result.shape[1]):
            temp = img_pad[i:i + m, j:j + n]
            img_result[i, j] = np.median(temp)
    return img_result

# 中值滤波器与高斯低通滤波器比较
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH03/Fig0335(a)(ckt_board_saltpep_prob_pt05).tif', 0)

kernel_19 = guass_kernel((19, 19), sigma=3)
kernel_19 = kernel_19 / kernel_19.sum()

img_19 = cv2.filter2D(img_ori, ddepth=-1, kernel=kernel_19)

img_median = median_filter(img_ori, (7, 7), mode='constant')

plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1), plt.imshow(img_ori, cmap='gray', vmin=0, vmax=255), plt.title(f'Original')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 2), plt.imshow(img_19, cmap='gray', vmin=0, vmax=255), plt.title(f'Guassian Filter 19x19')
plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 3), plt.imshow(img_median, cmap='gray', vmin=0, vmax=255), plt.title(f'Medial filter')
plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

jasneik

关注

1
点赞
踩
7

收藏

觉得还不错? 一键收藏
打赏
2
评论
第3章 Python 数字图像处理(DIP) - 灰度变换与空间滤波14 - 平滑低通滤波器 -高斯滤波器核的生成方法

目录平滑（低通）空间滤波器低通高斯滤波器核统计排序（非线性）滤波器平滑（低通）空间滤波器平滑（也称平均）空间滤波器用于降低灰度的急剧过渡在图像重取样之前平滑图像以减少混淆用于减少图像中无关细节平滑因灰度级数量不足导致的图像中的伪轮廓平滑核与一幅图像的卷积会模糊图像低通高斯滤波器核高斯核是唯一可分离的圆对称（也称各向同性，意味它们的响应与方向无关）核。w(s,t)=G(s,t)=Ke−s2+t22σ2(3.45)w(s,t) = G(s,t) = K e^{-\frac{s^2 + t^
复制链接

扫一扫