python图像处理《数字图像处理与python实现》读书笔记二:空间滤波

记录知识，形成体系

阿正的梦工坊

3054人浏览 · 2021-04-24 15:34:09

阿正的梦工坊 · 2021-04-24 15:34:09 发布

文章目录

重要！
第三章空间滤波

部分实验结果展示

彩色图像中值滤波实验结果

在这里插入图片描述
最大值滤波和最小值滤波实验结果

sobel边缘检测和梯度图像实验结果

灰度图像均值滤波实验结果
在这里插入图片描述

重要！

请读者写代码的时候，不时参考官方手册，可以节省很多时间。

scikit-image api参考手册：https://scikit-image.org/docs/0.18.x/api/api.html

matplotlib api 参考手册：https://matplotlib.org/stable/api/index.html

numpy api参考手册：https://www.numpy.org.cn/reference/routines/

scipy api参考手册：https://docs.scipy.org/doc/scipy/reference/api.html 本文主要用scipy库下的信号处理signal包。

写在前面：本文根据岳亚伟老师的《数字图像处理与python实现》一书整理而来。本文结合原书并动手进行实验，形成个人读书笔记，仅供个人学习使用，如有侵权，请联系删除。另，请读者尊重原作者知识成果。在此，对读者表示感谢。

第三章空间滤波

概览

空间滤波是在图像平面本身上逐像素地移动空间模板，同时空间模板与其覆盖的图像像素灰度值按预定义的关系进行运算。

模板也称为空间滤波器、核、掩模或窗口。

空间滤波一般用于去除图像噪声或增强图像细节，突出感兴趣信息，抑制无效信息，以改善人类的视觉效果或使图像更适合于特定的机器感知与分析。

空间滤波主要包括平滑处理和锐化处理两大类。

平滑处理主要用于去除图像中一些不重要的细节并减小噪声。

锐化处理主要用于突出图像中的细节，增强图像边缘。

为了达到较满意的图像增强效果，通常使用多种互补的空间滤波技术。

3.1 空间滤波基础

空间域指的是图像平面本身，是相对于变换域而言的。

空间域的图像处理是图像本身不进行频域变换，以图像中的像素为基础对图像进行处理。空间域的图像处理是在像素的邻域进行操作，如空间域平滑处理是通过像素的邻域平滑图像，空间域锐化处理是通过像素的邻域锐化图像。

频域的图像处理首先将图像变换到变换域，然后在频域进行处理，处理之后将其反变换至空间域。频域处理主要包括低通滤波和高通滤波。

低通滤波可以使低频信号正常通过，而高于所设定临界值的高频信号则被阻隔或减弱，可用于去除图像的噪声，相当于空间域的平滑处理。

高通滤波可以使高频信号正常通过，而低于所设定临界值的低频信号则被阻隔或减弱，可增强图像的边缘轮廓等高频信号，相当于空间域的锐化处理。

在频域处理中，滤波是指过滤一些频率分量，即通过一些频率分量，同时拒绝一些频率分量的通过。频域滤波器主要包括低通滤波器和高通滤波器。

滤波一词也可用于空间域，称为空间域滤波，即在空间域上直接对图像进行处理，实现类似于频域的平滑或锐化效果。

3.1.1 空间滤波的机理

空间滤波的机理就是在待处理图像上逐像素地移动模板，在每个像素点，滤波器的响应通过事先定义的关系计算。

若滤波器在图像像素上执行的是线性操作，则称为线性滤波器，否则称为非线性滤波器。

均值滤波器求解的是模板内像素灰度值的平均值，其是典型的线性滤波器。

统计排序滤波器是通过比较给定邻域内的灰度值大小实现的，原始数据与滤波结果是一种逻辑关系，如最大值滤波器、最小值滤波器、中值滤波器等，都是典型的非线性滤波器。

下面说明一下 $\times 3$ 模板的线性空间滤波器的机理，在图像的任一像素点 $(x, y)$ ,其灰度值为 $S (x, y)$ ，空间滤波器的响应 $T (x, y)$ 是模板系数与模板所覆盖的像素灰度值的乘积之和：

$T(x,y)=W(-1,-1)\times S(x-1,y-1)+W(-1,0)\times S(x-1,y)+...\\+W(0,0)\times S(x,y)+...+W(1,1)\times S(x+1,y+1)$

其中W表示3×3模板，S表示灰度图像，T表示线性空间滤波结果。显然，模板中心 $W (0, 0)$ 覆盖着像素点 $(x, y)$ .

在这里插入图片描述
来源：岳亚伟《数字图像处理与python实现》，下同

现将 $3\times 3$ 模板的线性空间滤波器推广到一般情形。对于一个大小为 $m\times n$ 的模板，由于主要关注奇数尺寸的模板，可以对m和n进行假定： $\times a + 1$ , $2\times b + 1$ ，其中a和b均为正整数。使用 $m\times n$ 的模板W对 $M\times N$ 大小的图像S进行线性空间滤波，得到滤波结果T图像。

在这里插入图片描述
其中 $(x, y)$ 表示图像中的任一像素点。

注意：式（3-1）是将模板与图像进行相关运算，而非卷积运算。本章的滤波器模板均是与图像进行相关运算，而非卷积运算。

实现空间滤波邻域处理时，需要考虑的一个问题是滤波中心靠近图像边界时如何计算空间滤波器的响应。考虑一个大小为m×n的模板，当模板中心距离图像左边界或右边界为(n-1)/2个像素时，该模板有一条边与图像左边界或右边界重合；当模板中心距离图像上边界或下边界为(m-1)/2个像素时，该模板有一条边与图像上边界或下边界重合。此时如果模板的中心继续向图像边界靠近，模板的行或列就会处于图像平面之外。

有多种方法可以解决上述问题，较简单的方法是将模板中心点的移动范围限制在距离图像左边界或右边界不小于(n-1)/2个像素，且距离图像上边界或下边界不小于(m-1)/2个像素。这种解决方法将使处理后的图像比原始图像稍小，可以将未被处理的像素点的灰度值直接复制到滤波结果处，以保持滤波结果与原始图像尺寸一致。

另一种解决方法是在图像的左边界和右边界以外各补上(n-1)/2列灰度为零的像素点（其灰度值可以为其他常量值，也可以是边界像素的灰度值），在图像的上边界和下边界以外各补上(m-1)/2行灰度为零的像素点，然后再进行滤波处理，处理后的图像与原始图像尺寸一致。

在这里插入图片描述
分析：滤波结果中，像素灰度值由模板系数与模板覆盖的像素灰度值的乘积之和求得。如滤波结果的第一个像素的灰度值2=1×0+0×0+0×0+0×0+0×1+0×2+0×0+0×0+2×1（相关运算）。即和滤波模板相同大小的矩阵，对应元素相乘并相加。

下面是对上述图片过程的代码实现：
linear_filtering_using_matrix.py

import numpy as np

def correl2d(img, window):
	m = window.shape[0] # 高度
	n = window.shape[1] # 宽度
	# 边界通过0灰度值进行填充
	img1 = np.zeros((img.shape[0] + m - 1, img.shape[1] + n -1))
	img1[(m-1) // 2 : (img.shape[0] + (m - 1) // 2), 
		(n - 1) // 2 : (img.shape[1] + (n - 1)//2)] = img
	img2 = np.zeros(img.shape)
	for i in range(img2.shape[0]):
		for j in range(img2.shape[1]):
			temp = img1[i : i + m, j : j + n] # 和模板大小相同的矩阵
			img2[i, j] = np.sum(np.multiply(temp, window)) # 矩阵对应元素相乘并求和
	return (img1, img2)

# window 表示滤波模板， img表示原始矩阵
window = np.array([[1, 0, 0], 
					[0, 0, 0], 
					[0, 0, 2]])
img = np.array([[1, 2, 1, 0, 2, 3],
				[0, 1, 1, 2, 0, 1],
				[3, 0, 2, 1, 2, 2], 
				[0, 1, 1, 0, 0, 1],
				[1, 1, 3, 2, 2, 0], 
				[0, 0, 1, 0, 1, 0]],)
# img1表示边界填充后的矩阵，img2表示空间滤波结果
img1, img2 = correl2d(img, window)
print("边界填充后的矩阵")
print(img1)
# print("\n")
print("空间滤波结果")
print(img2)

实验结果：

边界填充后的矩阵
[[0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 1. 2. 1. 0. 2. 3. 0.]
 [0. 0. 1. 1. 2. 0. 1. 0.]
 [0. 3. 0. 2. 1. 2. 2. 0.]
 [0. 0. 1. 1. 0. 0. 1. 0.]
 [0. 1. 1. 3. 2. 2. 0. 0.]
 [0. 0. 0. 1. 0. 1. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0.]]
空间滤波结果
[[2. 2. 4. 0. 2. 0.]
 [0. 5. 4. 5. 4. 2.]
 [2. 2. 1. 1. 4. 0.]
 [2. 9. 4. 6. 1. 2.]
 [0. 2. 1. 3. 0. 0.]
 [0. 1. 1. 3. 2. 2.]]

3.1.2 空间滤波器模板

若空间滤波器模板的系数从1开始进行索引，从左到右索引值递增，先索引第一行的每个模板系数，再依次索引下一行的每个模板系数，则3×3滤波模板的一种表示如图3-3所示。若用w向量表示滤波模板，z向量表示模板覆盖的像素的灰度值，则模板响应R可用式（3-2）表示。
在这里插入图片描述
其中， $w=[w_1,w_2, ..., w_{mn}]^T$ 是包含 $m\times n$ 个模板系数的向量形式的滤波器模板， $z=[z_1,z_2, ..., z_{mn}]^T$ 是滤波器模板覆盖的图像像素灰度值向量。

对于3×3模板，此时m=n=3，m×n=9，该模板响应为
在这里插入图片描述
其中w是包含9个模板系数的列向量，z是包含9个像素灰度值的列向量。

m×n的线性滤波器模板有m×n个模板系数，这m×n个系数的值决定了线性空间滤波器的功能。假设要实现3×3的平滑空间滤波器，较简单的方法是使得滤波器的系数均为1/9。3.2节中将会再次仔细讨论平滑空间滤波器，介绍多种平滑空间滤波模板。锐化空间滤波模板也将在3.3节进行详细讨论。

3.2 平滑处理

平滑处理常用于模糊处理和降低噪声。

平滑滤波器使用给定邻域内像素的平均灰度值或逻辑运算值代替原始图像中像素的灰度值，这种处理降低了图像灰度的“尖锐”变化。

然而，图像边缘也是由图像灰度尖锐变化带来的特性，因此平滑空间滤波器有边缘模糊化的负面效应。

平滑空间滤波器可分为平滑线性空间滤波器和平滑非线性空间滤波器。具有代表性的平滑非线性空间滤波器为统计排序滤波器。

3.2.1 平滑线性空间滤波器

平滑线性空间滤波器的输出是给定邻域内的像素灰度值的简单平均值或加权平均值。平滑线性空间滤波器有时也称为均值滤波器。

均值滤波器的一个重要应用是降低图像中的噪声。均值滤波器还有一个重要应用，去除图像的不相关细节，使不相关细节与背景糅合在一起，从而使感兴趣目标更加易于检测，此时模板的大小与不相关细节的尺寸有关。

在这里插入图片描述
上图展示了两个平滑线性空间滤波器模板。（a）是 $3\times 3$ 盒状滤波器模板；（b）是 $5\times 5$ 盒状滤波器模板.滤波器响应如下式：

R是由m×n大小的模板定义的均值滤波器的响应，该模板中的所有系数均为1/mn，这种滤波器也称为盒状滤波器，是最简单的均值滤波器。

盒状滤波是最简单的平滑空间滤波方法之一。使用3×3、5×5、9×9盒状滤波器对图像进行盒状滤波的代码如下。

import numpy as np
from scipy import signal
from skimage import data
from matplotlib import pyplot as plt

# 定义二维灰度图像的空间滤波函数
def correl2d(img, window):
	# 使用滤波器实现图像的空间相关
	# mode = 'same'表示输出尺寸等于输入尺寸
	# boundary = 'fill'表示滤波前，用常量值填充原始图像的边缘，默认常量值为0
	s = signal.correlate2d(img, window, mode='same', boundary='fill')
	return s.astype(np.uint8)
# img为原始图像
img = data.camera()
# 3*3盒状滤波模板
window_1 = np.ones((3, 3))/(3 ** 2)
# 5*5盒状滤波模板
window_2 = np.ones((5, 5))/(5 ** 2)
# 9*9盒状滤波模板
window_3 = np.ones((9, 9))/(9 ** 2)
# 生成滤波结果
new_img_1 = correl2d(img, window_1)
new_img_2 = correl2d(img, window_2)
new_img_3 = correl2d(img, window_3)
# 显示图像
plt.rcParams['font.sans-serif'] = ['SimHei'] # 中文
# 设置figsize，防止图片重叠
fig, axs = plt.subplots(2, 2,figsize=(6, 15))

axs[0, 0].imshow(img, cmap='gray')
axs[0, 0].set_title("摄影师原图")
axs[0, 1].imshow(new_img_1, cmap='gray')
axs[0, 1].set_title("3*3盒状滤波模板")
axs[1, 0].imshow(new_img_2, cmap='gray')
axs[1, 0].set_title("5*5盒状滤波模板")
axs[1, 1].imshow(new_img_3, cmap='gray')
axs[1, 1].set_title("9*9盒状滤波模板")
plt.show()

灰度图像均值滤波实验结果：
在这里插入图片描述

在这里插入图片描述

分析：
可以发现3×3盒状滤波结果比原始图像的湖面水波更加平滑。与原始图像相比，使用5×5盒状滤波模板所得的图像的湖面水波更加平滑，并且远处的风景更加模糊，同时图像中摄影师也被模糊了。与原始图像和其他滤波结果相比，9×9盒状滤波结果更模糊。使用盒状滤波器对该图像进行滤波，随着滤波模板的增大，在远处风景模糊化的同时，也将图像中的摄影师模糊化了。

下面对上述代码中的部分函数进行整理：
1） numpy.ones函数
回顾numpy ones函数
函数原型：

numpy.ones(shape, dtype=None, order=‘C’, *, like=None)[source]
Return a new array of given shape and type, filled with ones.

函数参数：

Parameters
shapeint or sequence of ints
Shape of the new array, e.g., (2, 3) or 2.

dtypedata-type, optional
The desired data-type for the array, e.g., numpy.int8. Default is numpy.float64.

order{‘C’, ‘F’}, optional, default: C
Whether to store multi-dimensional data in row-major (C-style) or column-major (Fortran-style) order in memory.

likearray_like
Reference object to allow the creation of arrays which are not NumPy arrays. 
If an array-like passed in as like supports the __array_function__ protocol, the result will be defined by it. 
In this case, it ensures the creation of an array object 
compatible with that passed in via this argument.

Note

The like keyword is an experimental feature pending on acceptance of NEP 35.

New in version 1.20.0.

Returns
outndarray
Array of ones with the given shape, dtype, and order.

基本用法：

np.ones(5)
array([1., 1., 1., 1., 1.])

np.ones((5,), dtype=int)
array([1, 1, 1, 1, 1])

np.ones((2, 1))
array([[1.],
       [1.]])
       
s = (2,2)
np.ones(s)
array([[1.,  1.],
       [1.,  1.]])

2）scipy.signal.correlate2d函数

函数原型：

scipy.signal.correlate2d(in1, in2, mode=‘full’, boundary=‘fill’, fillvalue=0)

主要功能：
求2个二维矩阵的互相关（交叉相关）

Cross-correlate two 2-dimensional arrays.
Cross correlate in1 and in2 with output size determined by mode, and boundary conditions determined by boundary and fillvalue.

函数参数：
在这里插入图片描述

来源：https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.correlate2d.html#scipy.signal.correlate2d

测试代码
test_correlate2d_function.py

Use 2D cross-correlation to find the location of a template in a noisy image 在有噪声的图像中找到模板的位置

from scipy import signal
from scipy import misc
import numpy as np
import matplotlib.pyplot as plt
from skimage import io

face = misc.face(gray=True) - misc.face(gray=True).mean()
template = np.copy(face[300:365, 670:750])  # right eye
template -= template.mean()
face = face + np.random.randn(*face.shape) * 50  # add noise
corr = signal.correlate2d(face, template, boundary='symm', mode='same')
y, x = np.unravel_index(np.argmax(corr), corr.shape)  # find the match


fig, (ax_orig, ax_template, ax_corr) = plt.subplots(3, 1,
                                                    figsize=(6, 15))
ax_orig.imshow(face, cmap='gray')
ax_orig.set_title('Original')
ax_orig.set_axis_off()
ax_template.imshow(template, cmap='gray')
ax_template.set_title('Template')
ax_template.set_axis_off()
ax_corr.imshow(corr, cmap='gray')
ax_corr.set_title('Cross-correlation')
ax_corr.set_axis_off()
ax_orig.plot(x, y, 'ro')
plt.show()

如果运行官网代码出现图片一闪而过的情况，请移步笔者的另一篇博文：python画图fig.show()一闪而过的解决方法
找到模板位置的测试结果
在这里插入图片描述

常用的均值滤波器是加权平均的，即在计算滤波器响应时邻域中某些像素的权重较大。

下图所示的加权平均滤波器模板，模板中心位置的系数最大，模板其他位置的系数与距离模板中心的距离成反比。用户可根据实际目标调整加权平均滤波器模板各系数的权重，加权平均滤波器模板比盒状滤波器模板更符合实际应用需求。

在这里插入图片描述

加权平均滤波器模板更合理，距离越近的像素点权重越大，距离越远的像素点权重越小。高斯分布显然是一种可取的权重分配模式。高斯分布是一种钟形曲线，越接近中心，取值越大，越远离中心，取值越小。

高斯滤波器是一类根据高斯函数的形状选择权值的线性平滑滤波器。高斯平滑滤波器对于抑制服从正态分布的噪声非常有效。

二维高斯函数如式（3-5）所示。
在这里插入图片描述

其中，σ是标准差。根据以上二维高斯函数可以生成高斯平滑滤波器模板。假设现在要生成高斯滤波器模板，且要求模板系数之和为1，则滤波器系数如式（3-6）所示。
在这里插入图片描述

模板大小为m×n，为了方便索引，使得m=2×a+1且n=2×b+1，其中a和b均为正整数。模板中心为W(0,0)，模板中心的系数值最大。假设现在要生成5×5的高斯平滑空间滤波模板，则需要计算W(-2,-2)、W(-2,-1)、…、W(0,0)、…、W(2,1)、W(2,2)的值。下图为标准差σ为1.0时的5×5高斯平滑滤波器模板。

在这里插入图片描述
高斯平滑滤波是一种应用较广泛的平滑空间滤波方法之一。使用3×3、5×5、9×9的高斯平滑滤波器对图像进行高斯平滑滤波的代码如下。

import numpy as np
from scipy import signal
from skimage import data
from matplotlib import pyplot as plt
import math

def correl2d(img, window):
	"""二维灰度图像的空间滤波函数"""
	# 使用滤波器实现图像的空间相关 
	# mode = 'same'表示输出尺寸等于输入尺寸	
	# boundary = 'fill' 表示滤波前，用常量值填充原始图像的边缘，默认常量值为0
	s = signal.correlate2d(img, window, mode='same', boundary='fill')
	return s.astype(np.uint8)

def gauss(i, j, sigma):
	"""定义二维高斯函数"""
	return 1 / (2 * math.pi * sigma ** 2) * math.exp(-(i ** 2 + j ** 2) 
		/ (2 * sigma ** 2))

def gauss_window(radius, sigma):
	"""定义radius * radius 的高斯平滑模板"""
	window = np.zeros((radius * 2 + 1, radius * 2 + 1))
	for i in range(-radius, radius + 1):
		for j in range(-radius, radius + 1):
			window[i + radius][j + radius] = gauss(i, j, sigma)
	return window / np.sum(window)
	
# img为原始图像
img = data.camera()
# 3*3 高斯平滑滤波模板
window_1 = gauss_window(3, 1.0)
# 5*5 高斯平滑滤波模板
window_2 = gauss_window(5, 1.0)
# 9*9 高斯平滑滤波模板
window_3 = gauss_window(9, 1.0)
# 生成滤波结果
new_img_1 = correl2d(img, window_1)
new_img_2 = correl2d(img, window_2)
new_img_3 = correl2d(img, window_3)
# 显示图像
plt.rcParams['font.sans-serif'] = ['SimHei'] # 中文
fig, axs = plt.subplots(2, 2)
axs[0, 0].imshow(img,)
axs[0, 0].set_title("摄影师原图")
axs[0, 1].imshow(new_img_1,)
axs[0, 1].set_title("3*3 高斯平滑滤波模板")
axs[1, 0].imshow(new_img_2, )
axs[1, 0].set_title("5*5 高斯平滑滤波模板")
axs[1, 1].imshow(new_img_3,)
axs[1, 1].set_title("9*9 高斯平滑滤波模板")
plt.tight_layout()
plt.show()

灰度图像高斯滤波实验结果

在这里插入图片描述

对原始图像分别进行3×3高斯平滑滤波、5×5高斯平滑滤波、9×9高斯平滑滤波的结果如图所示，可以发现随着高斯滤波模板的增大，滤波结果越来越平滑。

对比均值滤波和高斯滤波，可以发现使用相同尺寸的模板，高斯滤波后图像被平滑的程度较低。高斯滤波的输出是邻域像素的加权平均，同时距离中心越近的像素权重越大。因此，与盒状滤波相比，高斯滤波的平滑效果更柔和，图像中感兴趣目标的细节保留得更好。高斯滤波效果与标准差、模板尺寸有关，可以调节标准差和模板尺寸，观察高斯平滑滤波的效果。

3.2.2 统计排序滤波器

统计排序滤波器是典型的非线性平滑滤波器，首先对模板覆盖的像素的灰度值进行排序，选择有代表性的灰度值作为统计排序滤波器的响应。

典型的统计排序滤波器包括最大值滤波器、中值滤波器和最小值滤波器。

最大值滤波器是用像素邻域内的最大值代替该像素的灰度值，主要用于寻找最亮点。
中值滤波器是用像素邻域内的中间值代替该像素的灰度值，主要用于降噪。
最小值滤波器是用像素邻域内的最小值代替该像素的灰度值，主要用于寻找最暗点。

在统计排序滤波器中，中值滤波器的应用最广。

对于一定类型的随机噪声，中值滤波器的降噪效果较好，同时比相同尺寸的均值滤波器模糊程度明显要低。中值滤波器对处理脉冲噪声（也称椒盐噪声）非常有效，因为中值滤波器取中值作为滤波结果，可以很好地去除滤波器覆盖的邻域中的一些黑点或者白点。

中值滤波器首先对模板覆盖的像素邻域内的所有灰度值进行排序，找到邻域的中间值，用这个中间值作为中值滤波器的响应。假设3×3中值滤波模板覆盖的像素灰度值为（2、3、0、10、9、1、7、5、3），排序结果为（0、1、2、3、3、5、7、9、10），中间值为3，则该邻域的中值滤波结果为3。中值滤波器使得图像中突出的亮点（暗点）更像它周围的值，以消除孤立的亮点（暗点），从而实现对图像的平滑。

为了观察中值滤波的降噪效果，首先对宇航员的灰度图像加入椒盐噪声，然后使用3×3的中值滤波器对图像进行中值滤波。中值滤波代码如下。

from scipy import ndimage
from skimage import data, util
from matplotlib import pyplot as plt

# img为原始图像
img = data.astronaut()[:, :, 0]
# 对图像加上椒盐噪声
noise_img = util.random_noise(img, mode='s&p', seed=None, clip=True)
# 中值滤波
n = 3
new_img = ndimage.median_filter(noise_img, (n, n))
# 显示图像
plt.rcParams['font.sans-serif'] = ['SimHei'] # 中文
fig, axs = plt.subplots(1, 3)
axs[0].imshow(img, cmap='gray')
axs[0].set_title("宇航员原图")
axs[1].imshow(noise_img, cmap='gray')
axs[1].set_title("加噪结果")
axs[2].imshow(new_img, cmap='gray')
axs[2].set_title("降噪结果")

plt.tight_layout()
plt.show()

灰度图像中值滤波实验结果

在这里插入图片描述

观察上图，可以发现中值滤波后的图像与原始图像非常接近。中值滤波可以很好地去除随机椒盐噪声。中值滤波的应用较广泛。在实践中常使用中值滤波器对图像进行降噪处理。

以上空间滤波的例子都是在灰度图像上进行的，其实也可以对彩色图像进行空间滤波。对RGB图像的空间滤波相当于分别对R、G、B三通道的图像进行空间滤波。对宇航员的彩色图像加入椒盐噪声，然后使用3×3的中值滤波器对彩色图像进行平滑空间滤波。彩色图像中值滤波代码如下。

import numpy as np
from scipy import ndimage
from skimage import data, util
from matplotlib import pyplot as plt

# img为原始图像
img = data.astronaut()
noise_img = np.zeros(img.shape)
new_img = np.zeros(img.shape)

for i in range(3):
	gray_img = data.astronaut()[:, :, i]
	# 对图像加上椒盐噪声
	noise_img[:, :, i] = util.random_noise(gray_img, mode='s&p', 
											seed=None, clip=True)
	n = 3
	new_img[:, :, i] = ndimage.median_filter(noise_img[:, :, i], (n, n))

# 显示图像
plt.rcParams['font.sans-serif'] = ['SimHei'] # 中文
fig, axs = plt.subplots(1, 3)
axs[0].imshow(img, cmap='gray')
axs[0].set_title("宇航员原图")
axs[1].imshow(noise_img, cmap='gray')
axs[1].set_title("加噪结果")
axs[2].imshow(new_img, cmap='gray')
axs[2].set_title("降噪结果")

plt.tight_layout()
plt.show()

彩色图像中值滤波实验结果

在这里插入图片描述

大图显示
在这里插入图片描述

最大值滤波器是将邻域内的像素灰度值进行从小到大的排序，用序列的最后一个值（即最大值）代替该像素的灰度值，对于发现图像最亮点非常有效，可有效降低胡椒噪声。最小值滤波器用序列的最小值代替该像素的灰度值，对于发现图像最暗点非常有效，可有效降低盐粒噪声。

对图像加入噪声，然后分别使用3×3的最大值滤波器和最小值滤波器对图像进行空间滤波，以观察最大值滤波和最小值滤波效果。最大值滤波和最小值滤波代码如下。

from scipy import ndimage
from skimage import data, util
from matplotlib import pyplot as plt
import matplotlib.gridspec as gridspec

img = data.astronaut()[:, :, 0]
# 对图像加入胡椒噪声
pepper_img = util.random_noise(img, mode='pepper', seed=None, clip=True)
# 对图像加入盐粒噪声
salt_img = util.random_noise(img, mode='salt', seed=None, clip=True)
n = 3
# 最大值滤波
max_img = ndimage.maximum_filter(pepper_img, (n, n))
# 最小值滤波
min_img = ndimage.minimum_filter(salt_img, (n, n))
# 显示图像
plt.rcParams['font.sans-serif'] = ['SimHei'] # 中文
gs = gridspec.GridSpec(2, 6)
gs.update(wspace=0.8)
ax1 = plt.subplot(gs[0, :2])
ax1.imshow(img,cmap='gray')
ax1.set_title("宇航员原图")

ax2 = plt.subplot(gs[0, 2:4])
ax2.imshow(pepper_img,cmap='gray')
ax2.set_title("加胡椒噪声图像")


ax3 = plt.subplot(gs[0, 4:6])
ax3.imshow(salt_img,cmap='gray')
ax3.set_title("加盐粒噪声图像")

ax4 = plt.subplot(gs[1, 1:3])
ax4.imshow(max_img,cmap='gray')
ax4.set_title("最大值滤波结果")

ax5 = plt.subplot(gs[1, 3:5])
ax5.imshow(min_img,cmap='gray')
ax5.set_title("最小值滤波结果")

plt.show()

最大值滤波和最小值滤波实验结果：
在这里插入图片描述

最大值滤波和最小值滤波分析：

可以发现最大值滤波结果比原始图像亮，最小值滤波结果比原始图像暗。最大值滤波对于去除胡椒噪声非常有效，最小值滤波对于去除盐粒噪声非常有效。

其实，这里笔者做实验时遇到问题，然后得以解决：python绘图subplot绘制5幅图：以2行绘制，首行3幅图，次行2幅图居中(内含绘制3幅图简单版)

下面给出测试绘制3幅图的代码：

import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
gs = gridspec.GridSpec(2, 4)
gs.update(wspace=0.5)
ax1 = plt.subplot(gs[0, :2], )
ax2 = plt.subplot(gs[0, 2:])
ax3 = plt.subplot(gs[1, 1:3])
plt.show()

绘制三幅图结果：
在这里插入图片描述
接着是subplot绘制五幅图的代码

import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
gs = gridspec.GridSpec(2, 6)
gs.update(wspace=0.8)
ax1 = plt.subplot(gs[0, :2])
ax2 = plt.subplot(gs[0, 2:4])
ax3 = plt.subplot(gs[0, 4:6])
ax4 = plt.subplot(gs[1, 1:3])
ax5 = plt.subplot(gs[1, 3:5])
plt.show()

绘制5幅图的结果
在这里插入图片描述

平滑处理是基本的图像处理方法之一。在实践中，盒状滤波器、高斯平滑滤波器、中值滤波器等都是常用的滤波器，常用于图像降噪或者图像预处理。在实践中，对图像进行平滑处理时，选用何种滤波器以及滤波模板的大小须结合实际目标。

3.3 锐化处理

锐化处理的目的是增强图像中目标的细节、边缘、轮廓和其他灰度突变，削弱了灰度变化缓慢的区域。

由于微分是对函数的局部变化率的一种描述，因此图像锐化算法的实现可基于空间微分。

图像平滑处理有边缘和细节模糊的负面效应。图像平滑和图像锐化在逻辑上是相反的操作，因此也可以使用原始图像减去平滑处理后的图像实现锐化处理，这称为反锐化掩蔽。

3.3.1 一阶微分算子

对任意一阶微分的定义都必须满足以下两点：在灰度不变的区域微分值为0；在灰度变化的区域微分值非0。由于处理的是离散情况，微分用差分近似。对于一维函数f(x)，其一阶微分的基本定义是：
在这里插入图片描述

一维微分可以用导数符号，这里使用偏导数符号是为了方便扩展至二维微分。二维图像f(x,y)将沿着两个空间坐标轴求解一阶微分，即分别求解f(x,y)对x和y的偏导数：
在这里插入图片描述
二维图像的梯度是一个二维列向量，可表示为：

梯度的幅值M(x,y)可表示为：

M(x,y)是梯度向量方向变化率在(x,y)处的值。梯度向量的两个分量都是一阶偏导数，说明两个分量都是线性算子。梯度向量的幅值进行了求平方、求平方根操作，是非线性算子。

为了方便描述，将像素点 $(x, y)$ 的3×3邻域的灰度值表示为 $z_1,z_2,..,z_9$ ,其中 $z_5$ 表示位于 $(x, y)$ 的像素的灰度值 $f (x, y)$ ， $z_1$ 表示 $f (x - 1, y - 1)$ ， $z_6$ 表示 $f (x + 1, y)$ ， $z_8$ 表示 $f (x, y + 1)$ ， $z_9$ 表示 $f (x + 1, y + 1)$ ，图像的3×3邻域如下图所示。

在这里插入图片描述

二维图像的最简单的一阶微分近似是： $g_x=z_8-z_5$ 和 $g_y=z_6-z_5$ 。

根据罗伯特（Roberts）的观点，边缘探测器应具有以下特性：产生的边缘应清晰；背景应尽可能减少噪声；边缘强度应尽可能接近人类的感知。考虑到图像边界的拓扑结构性，罗伯特提出两个交叉差分表示 $g_x$ 和 $g_y$ :

$g_x=z_9-z_5,g_y=z_8-z_6$

罗伯特提出两个交叉差分也称为罗伯特交叉梯度算子，如下图所示。罗伯特正对角线梯度算子和负对角线梯度算子分别如下图（a）和（b）所示。梯度幅度M(x,y)为：
在这里插入图片描述

在这里插入图片描述

求罗伯特边缘图像和梯度图像的代码如下。该代码调用了skimage.filters中的roberts_pos_diag()、roberts_neg_diag()、roberts()方法，分别求取罗伯特正对角线边缘图像、罗伯特负对角线边缘图像、罗伯特梯度图像。

roberts边缘检测和梯度图像实验结果：
在这里插入图片描述
将罗伯特正对角线梯度算子作用于原始图像，可得到罗伯特正对角线边缘图像，将罗伯特负对角线梯度算子作用于原始图像，可得到罗伯特负对角线边缘图像，roberts梯度图像可根据梯度幅度公式（平方和的平方根）求得。

roberts边缘检测和梯度图像源代码

from skimage import data, filters
from matplotlib import pyplot as plt

img = data.camera()
# robert梯度交叉算子
img_robert_pos = filters.roberts_pos_diag(img)
img_robert_neg = filters.roberts_neg_diag(img)
img_robert = filters.roberts(img)

# 显示图像
plt.rcParams['font.sans-serif'] = ['SimHei'] # 中文
fig, axs = plt.subplots(2, 2)
axs[0, 0].imshow(img, cmap='gray')
axs[0, 0].set_title("摄影师原图")
axs[0, 1].imshow(img_robert_pos, cmap='gray')
axs[0, 1].set_title("roberts正对角线边缘图像")
axs[1, 0].imshow(img_robert_neg, cmap='gray')
axs[1, 0].set_title("roberts负对角线边缘图像")
axs[1, 1].imshow(img_robert, cmap='gray')
axs[1, 1].set_title("roberts梯度图像")
plt.tight_layout()
plt.show()

关于filters.roberts函数的官方文档：https://scikit-image.org/docs/stable/api/skimage.filters.html#skimage.filters.roberts

函数详解如下：
在这里插入图片描述

roberts函数源代码：https://github.com/scikit-image/scikit-image/blob/v0.18.0/skimage/filters/edges.py#L552-L585

def roberts(image, mask=None):
    """Find the edge magnitude using Roberts' cross operator.

    Parameters
    ----------
    image : 2-D array
        Image to process.
    mask : 2-D array, optional
        An optional mask to limit the application to a certain area.
        Note that pixels surrounding masked regions are also masked to
        prevent masked regions from affecting the result.

    Returns
    -------
    output : 2-D array
        The Roberts' Cross edge map.

    See also
    --------
    sobel, scharr, prewitt, feature.canny

    Examples
    --------
    >>> from skimage import data
    >>> camera = data.camera()
    >>> from skimage import filters
    >>> edges = filters.roberts(camera)

    """
    check_nD(image, 2)
    out = np.sqrt(roberts_pos_diag(image, mask) ** 2 +
                  roberts_neg_diag(image, mask) ** 2)
    out /= np.sqrt(2)
    return out

可以发现，roberts函数就是对梯度幅度的计算。

roberts_pos_diag函数源代码如下：

用Roberts交叉算子求图像的交叉边。将核函数应用于输入图像，产生一个方向上的梯度分量的单独测量值。

def roberts_pos_diag(image, mask=None):
    """Find the cross edges of an image using Roberts' cross operator.
    The kernel is applied to the input image to produce separate measurements
    of the gradient component one orientation.
    Parameters
    ----------
    image : 2-D array
        Image to process.
    mask : 2-D array, optional
        An optional mask to limit the application to a certain area.
        Note that pixels surrounding masked regions are also masked to
        prevent masked regions from affecting the result.
    Returns
    -------
    output : 2-D array
        The Robert's edge map.
    Notes
    -----
    We use the following kernel::
      1   0
      0  -1
    """
    check_nD(image, 2)
    image = img_as_float(image)
    result = convolve(image, ROBERTS_PD_WEIGHTS)
    return _mask_filter_result(result, mask)

用到的ROBERTS_PD_WEIGHTS定义如下：

# 2D-only filter weights
ROBERTS_PD_WEIGHTS = np.array([[1, 0],
                               [0, -1]], dtype=np.double)
ROBERTS_ND_WEIGHTS = np.array([[0, 1],
                               [-1, 0]], dtype=np.double)

用到的check_nD函数源代码：https://github.com/scikit-image/scikit-image/blob/main/skimage/_shared/utils.py

def check_nD(array, ndim, arg_name='image'):
    """
    Verify an array meets the desired ndims and array isn't empty.
    Parameters
    ----------
    array : array-like
        Input array to be validated
    ndim : int or iterable of ints
        Allowable ndim or ndims for the array.
    arg_name : str, optional
        The name of the array in the original function.
    """
    array = np.asanyarray(array)
    msg_incorrect_dim = "The parameter `%s` must be a %s-dimensional array"
    msg_empty_array = "The parameter `%s` cannot be an empty array"
    if isinstance(ndim, int):
        ndim = [ndim]
    if array.size == 0:
        raise ValueError(msg_empty_array % (arg_name))
    if not array.ndim in ndim:
        raise ValueError(msg_incorrect_dim % (arg_name, '-or-'.join([str(n) for n in ndim])))

当然，上面调用了很多函数，读者感兴趣可以参考上述链接所指向的源码仓库。

使用罗伯特交叉梯度算子可得到梯度图像 $M (x, y)$ ，那如何对图像进行锐化，以增强图像的边缘呢？ 将梯度图像以一定比例叠加到原始图像 $f (x, y)$ ，即可得到锐化图像，如式（3-7）所示:

在这里插入图片描述
其中c为锐化强度系数.

由于奇数模板有对称中心，更易于实现，因此一般更注重奇数模板。罗伯特交叉梯度算子是2×2偶数模板，而我们更关注3×3奇数模板。使用上文 $z_1,...,z_9$ 表示法所用的3×3邻域表示方法，可以对3×3模板的 $g_x$ 和 $g_y$ 进行近似表达:
在这里插入图片描述

式（3-8）可用下图（a）所示的水平索贝尔（Sobel）算子和图（b）所示的竖直索贝尔算子实现。
在这里插入图片描述

sobel边缘检测和梯度图像实验结果
在这里插入图片描述

sobel边缘检测和梯度图像实验分析
将水平索贝尔算子作用于原始图像，可得到水平索贝尔边缘图像，其水平边缘较为明显。将竖直索贝尔算子作用于原始图像，可得到竖直索贝尔边缘图像，其竖直边缘较为明显。在以上代码中，索贝尔梯度图像的获取调用了skimage.filters中的索贝尔方法，索贝尔梯度图像也可以由第二幅图和第三幅图的平方和的平方根求得。

sobel边缘检测和梯度图像实验源代码

from skimage import data, filters
from matplotlib import pyplot as plt

img = data.camera()
# sobel 算子
img_sobel_h = filters.sobel_h(img)
img_sobel_v = filters.sobel_v(img)
img_sobel = filters.sobel(img)

# 显示图像
plt.rcParams['font.sans-serif'] = ['SimHei'] # 中文
fig, axs = plt.subplots(2, 2)
axs[0, 0].imshow(img, cmap='gray')
axs[0, 0].set_title("摄影师原图")
axs[0, 1].imshow(img_sobel_h, cmap='gray')
axs[0, 1].set_title("水平sobel边缘图像")
axs[1, 0].imshow(img_sobel_v, cmap='gray')
axs[1, 0].set_title("竖直sobel边缘图像")
axs[1, 1].imshow(img_sobel, cmap='gray')
axs[1, 1].set_title("sobel梯度图像")
plt.tight_layout()
plt.show()

关于sobel算子的函数

在这里插入图片描述

来源：官方文档

sobel函数 api ：https://scikit-image.org/docs/0.18.x/api/skimage.filters.html#skimage.filters.sobel
sobel函数源代码：
github仓库 https://github.com/scikit-image/scikit-image/blob/v0.18.0/skimage/filters/edges.py#L244-L271

def sobel(image, mask=None, *, axis=None, mode='reflect', cval=0.0):
    """Find edges in an image using the Sobel filter.
    Parameters
    ----------
    image : array
        The input image.
    mask : array of bool, optional
        Clip the output image to this mask. (Values where mask=0 will be set
        to 0.)
    axis : int or sequence of int, optional
        Compute the edge filter along this axis. If not provided, the edge
        magnitude is computed. This is defined as::
            sobel_mag = np.sqrt(sum([sobel(image, axis=i)**2
                                     for i in range(image.ndim)]) / image.ndim)
        The magnitude is also computed if axis is a sequence.
    mode : str or sequence of str, optional
        The boundary mode for the convolution. See `scipy.ndimage.convolve`
        for a description of the modes. This can be either a single boundary
        mode or one boundary mode per axis.
    cval : float, optional
        When `mode` is ``'constant'``, this is the constant used in values
        outside the boundary of the image data.
    Returns
    -------
    output : array of float
        The Sobel edge map.
    See also
    --------
    scharr, prewitt, canny
    References
    ----------
    .. [1] D. Kroon, 2009, Short Paper University Twente, Numerical
           Optimization of Kernel Based Image Derivatives.
    .. [2] https://en.wikipedia.org/wiki/Sobel_operator
    Examples
    --------
    >>> from skimage import data
    >>> from skimage import filters
    >>> camera = data.camera()
    >>> edges = filters.sobel(camera)
    """
    image = img_as_float(image)
    output = _generic_edge_filter(image, smooth_weights=SOBEL_SMOOTH,
                                  axis=axis, mode=mode, cval=cval)
    output = _mask_filter_result(output, mask)
    return output

用到的常量SOBEL_SMOOTH：

# n-dimensional filter weights
SOBEL_EDGE = np.array([1, 0, -1])
SOBEL_SMOOTH = np.array([1, 2, 1]) / 4
HSOBEL_WEIGHTS = SOBEL_EDGE.reshape((3, 1)) * SOBEL_SMOOTH.reshape((1, 3))
VSOBEL_WEIGHTS = HSOBEL_WEIGHTS.T

补充numpy.array的用法

>>>np.array([1, 2, 3])
array([1, 2, 3])

其中用到的_generic_edge_filter函数源代码：
该函数用于：通过沿一个维度应用边缘权重和沿所有其他维度应用平滑权重来生成滤波器。如果没有给定轴，或者给定了一个轴元组，则沿所有轴依次计算滤波器，并将幅值计算为所有轴的平均幅值平方根。


def _generic_edge_filter(image, *, smooth_weights, edge_weights=[1, 0, -1],
                         axis=None, mode='reflect', cval=0.0, mask=None):
    """Apply a generic, n-dimensional edge filter.
    The filter is computed by applying the edge weights along one dimension
    and the smoothing weights along all other dimensions. If no axis is given,
    or a tuple of axes is given the filter is computed along all axes in turn,
    and the magnitude is computed as the square root of the average square
    magnitude of all the axes.
    Parameters
    ----------
    image : array
        The input image.
    smooth_weights : array of float
        The smoothing weights for the filter. These are applied to dimensions
        orthogonal to the edge axis.
    edge_weights : 1D array of float, optional
        The weights to compute the edge along the chosen axes.
    axis : int or sequence of int, optional
        Compute the edge filter along this axis. If not provided, the edge
        magnitude is computed. This is defined as::
            edge_mag = np.sqrt(sum([_generic_edge_filter(image, ..., axis=i)**2
                                    for i in range(image.ndim)]) / image.ndim)
        The magnitude is also computed if axis is a sequence.
    mode : str or sequence of str, optional
        The boundary mode for the convolution. See `scipy.ndimage.convolve`
        for a description of the modes. This can be either a single boundary
        mode or one boundary mode per axis.
    cval : float, optional
        When `mode` is ``'constant'``, this is the constant used in values
        outside the boundary of the image data.
    """
    ndim = image.ndim
    if axis is None:
        axes = list(range(ndim))
    elif np.isscalar(axis):
        axes = [axis]
    else:
        axes = axis
    return_magnitude = (len(axes) > 1)

    output = np.zeros(image.shape, dtype=float)

    for edge_dim in axes:
        kernel = _reshape_nd(edge_weights, ndim, edge_dim)
        smooth_axes = list(set(range(ndim)) - {edge_dim})
        for smooth_dim in smooth_axes:
            kernel = kernel * _reshape_nd(smooth_weights, ndim, smooth_dim)
        ax_output = ndi.convolve(image, kernel, mode=mode)
        if return_magnitude:
            ax_output *= ax_output
        output += ax_output

    if return_magnitude:
        output = np.sqrt(output) / np.sqrt(ndim)
    return output

而后的sobel_h 函数源码

def sobel_h(image, mask=None):
    """Find the horizontal edges of an image using the Sobel transform.

    Parameters
    ----------
    image : 2-D array
        Image to process.
    mask : 2-D array, optional
        An optional mask to limit the application to a certain area.
        Note that pixels surrounding masked regions are also masked to
        prevent masked regions from affecting the result.

    Returns
    -------
    output : 2-D array
        The Sobel edge map.

    Notes
    -----
    We use the following kernel::

      1   2   1
      0   0   0
     -1  -2  -1

    """
    check_nD(image, 2)
    return sobel(image, mask=mask, axis=0)

咱们代码的中的调用：

img_sobel_h = filters.sobel_h(img)

当然，sobel_v用法雷同。

3.3.2 二阶微分算子

任意二阶微分的定义都必须满足以下3点：

在灰度不变的区域微分值为0；
在灰度台阶或斜坡的起点处微分值非0；
沿着斜坡的微分值为0。

由于我们处理的是离散情况，因此微分用差分近似。对于一维函数f(x)，其二阶微分的基本定义是：
在这里插入图片描述

对于二维图像f(x,y)，将沿着两个空间坐标轴求解二阶微分：
在这里插入图片描述

则拉普拉斯算子可表示为：

在这里插入图片描述

拉普拉斯算子可用下图（a）实现。上述拉普拉斯变换未考虑对角线元素，可以对其添加对角线元素项，并且更改中心项的系数，以保证模板系数和为0，从而保证灰度恒定区域的微分值为0。扩展的拉普拉斯算子如式（3-9）所示：
在这里插入图片描述
扩展的拉普拉斯算子如图（b）所示。图（a）和图（b）的中心系数均为负数，在实践中并不经常使用，实践中常使用的两个拉普拉斯算子是图（a）和图（b）分别乘以系数-1所得，如图（c）和（d）所示。

综上：Laplace算子常用的是中心系数为正数的版本！！

0 -1 0
-1 4 -1
0 -1 0

在这里插入图片描述

对原始图像使用拉普拉斯算子进行空间滤波可得到拉普拉斯图像，将拉普拉斯图像以一定比例叠加到原始图像（该比例系数的符号与拉普拉斯模板中心系数的符号一致），可对原始图像进行拉普拉斯锐化增强。拉普拉斯锐化增强的代码如下。

二阶微分算子laplace图像增强实验结果
在这里插入图片描述

补充：一阶微分算子sobel边缘检测和梯度图像实验结果
在这里插入图片描述
laplace图像增强结果分析
对原始图像使用拉普拉斯算子进行空间滤波可得到拉普拉斯图像。将原始图像和拉普拉斯图像以一定比例叠加，可得到拉普拉斯锐化增强图像。对比，可发现一阶微分更加突出图像的边缘，二阶微分对灰度变化强烈的地方更敏感，更加突出图像的纹理结构。

laplace图像增强源代码

from skimage import data, filters
from matplotlib import pyplot as plt

img = data.camera()
# laplace 算子
img_laplace = filters.laplace(img, ksize=3, mask=None)
img_enhance = img + img_laplace

# 显示图像
plt.rcParams['font.sans-serif'] = ['SimHei'] # 中文
fig, axs = plt.subplots(1, 3)
axs[0].imshow(img, cmap='gray')
axs[0].set_title("摄影师原图")
axs[1].imshow(img_laplace, cmap='gray')
axs[1].set_title("laplace图像")
axs[2].imshow(img_enhance, cmap='gray')
axs[2].set_title("laplace锐化增强图像")

plt.tight_layout()
plt.show()

对于filters.laplace函数的详解
在这里插入图片描述

laplace函数源码

def laplace(image, ksize=3, mask=None):
    """Find the edges of an image using the Laplace operator.

    Parameters
    ----------
    image : ndarray
        Image to process.
    ksize : int, optional
        Define the size of the discrete Laplacian operator such that it
        will have a size of (ksize,) * image.ndim.
    mask : ndarray, optional
        An optional mask to limit the application to a certain area.
        Note that pixels surrounding masked regions are also masked to
        prevent masked regions from affecting the result.

    Returns
    -------
    output : ndarray
        The Laplace edge map.

    Notes
    -----
    The Laplacian operator is generated using the function
    skimage.restoration.uft.laplacian().

    """
    image = img_as_float(image)
    # Create the discrete Laplacian operator - We keep only the real part of
    # the filter
    _, laplace_op = laplacian(image.ndim, (ksize,) * image.ndim)
    result = convolve(image, laplace_op)
    return _mask_filter_result(result, mask)

其中laplacian函数的源代码：
https://github.com/scikit-image/scikit-image/blob/main/skimage/restoration/uft.py

函数功能：返回拉普拉斯函数的传递函数。

def laplacian(ndim, shape, is_real=True):
    """Return the transfer function of the Laplacian.
    Laplacian is the second order difference, on row and column.
    Parameters
    ----------
    ndim : int
        The dimension of the Laplacian.
    shape : tuple
        The support on which to compute the transfer function.
    is_real : boolean, optional
       If True (default), imp_resp is assumed to be real-valued and
       the Hermitian property is used with rfftn Fourier transform
       to return the transfer function.
    Returns
    -------
    tf : array_like, complex
        The transfer function.
    impr : array_like, real
        The Laplacian.
    Examples
    --------
    >>> tf, ir = laplacian(2, (32, 32))
    >>> np.all(ir == np.array([[0, -1, 0], [-1, 4, -1], [0, -1, 0]]))
    True
    >>> np.all(tf == ir2tf(ir, (32, 32)))
    True
    """
    impr = np.zeros([3] * ndim)
    for dim in range(ndim):
        idx = tuple([slice(1, 2)] * dim +
                    [slice(None)] +
                    [slice(1, 2)] * (ndim - dim - 1))
        impr[idx] = np.array([-1.0,
                              0.0,
                              -1.0]).reshape([-1 if i == dim else 1
                                              for i in range(ndim)])
    impr[(slice(1, 2), ) * ndim] = 2.0 * ndim
    return ir2tf(impr, shape, is_real=is_real), impr

3.3.3 反锐化掩蔽

图像平滑处理有边缘和细节模糊的负面效应，因此可用原始图像减去平滑处理后的图像实现锐化处理，称为反锐化掩蔽。

反锐化掩蔽处理包括3个步骤：

通过平滑滤波得到模糊图像;
从原始图像中减去模糊图像得到差值图像;
将差值图像叠加到原始图像中。

原始图像 $f (x, y)$ 平滑处理所得的模糊图像为 $s (x, y)$ ，用原始图像减去模糊图像得到差值图像 $d (x, y)$ ，如式（3-10）所示。

在这里插入图片描述

最后，将差值图像以一定比例叠加到原始图像，如式（3-11）所示:
在这里插入图片描述
其中，权重系数为c(c≥0)。c=1时称为反锐化掩蔽，c>1时称为高提升滤波，c<1时不强调反锐化掩蔽效果。

反锐化掩蔽实验

反锐化掩蔽实验结果
在这里插入图片描述

原始图像进行3×3盒状滤波产生了模糊图像，原始图像与模糊图像做差值运算，得到差值图像。差值图像中的边缘信息较为丰富。原始图像与差值图像相叠加可得到锐化增强图像.读者可以修改以上代码中原始图像与差值图像的叠加比例，以观察图像反锐化掩蔽和高提升滤波的效果。

反锐化掩蔽实验代码

import numpy as np
from scipy import signal
from skimage import data
from matplotlib import pyplot as plt
def correl2d(img, window):
	"""
	定义二维灰度图像的空间滤波函数
	使用滤波器实现图像的空间相关
	mode = 'same'表示输出尺寸等于输入尺寸
	boundary = 'fill'表示滤波前，用常量值填充原始图像的边缘，默认常量值为0
	"""
	s = signal.correlate2d(img, window, mode='same', boundary='fill')
	return s

img = data.camera()
# 3*3盒状滤波模板
window = np.ones((3, 3)) / (3 ** 2)
img_blur = correl2d(img, window)
img_edge = img - img_blur
img_enhance = img + img_edge

# 显示图像
plt.rcParams['font.sans-serif'] = ['SimHei'] # 中文
fig, axs = plt.subplots(2, 2)
axs[0, 0].imshow(img, cmap='gray')
axs[0, 0].set_title("摄影师原图")
axs[0, 1].imshow(img_blur, cmap='gray')
axs[0, 1].set_title("模糊图像")
axs[1, 0].imshow(img_edge, cmap='gray')
axs[1, 0].set_title("差值图像")
axs[1, 1].imshow(img_enhance,cmap='gray')
axs[1, 1].set_title("锐化增强图像")
plt.tight_layout()
plt.show()

可调整权重系数c的代码

import numpy as np
from scipy import signal
from skimage import data
from matplotlib import pyplot as plt
PICTURE_ENHANCE_RATIO= 1 # 权重系数：为1表示反锐化掩蔽
def correl2d(img, window):
	s = signal.correlate2d(img, window, mode='same', boundary='fill')
	return s

img = data.camera()
# 3*3盒状滤波模板
window = np.ones((3, 3)) / (3 ** 2)
img_blur = correl2d(img, window)
img_edge = img - img_blur
img_enhance = img + PICTURE_ENHANCE_RATIO * img_edge

# 显示图像
plt.rcParams['font.sans-serif'] = ['SimHei'] # 中文
fig, axs = plt.subplots(2, 2)
axs[0, 0].imshow(img, cmap='gray')
axs[0, 0].set_title("摄影师原图")
axs[0, 1].imshow(img_blur, cmap='gray')
axs[0, 1].set_title("模糊图像")
axs[1, 0].imshow(img_edge, cmap='gray')
axs[1, 0].set_title("差值图像")
axs[1, 1].imshow(img_enhance,cmap='gray')
axs[1, 1].set_title("锐化增强图像")
plt.tight_layout()
plt.show()

c = 2时实验结果
在这里插入图片描述
c = 20时实验结果

在这里插入图片描述
c = 0.5时实验结果

上述c的含义：权重系数为c(c≥0)。c=1时称为反锐化掩蔽，c>1时称为高提升滤波，c<1时不强调反锐化掩蔽效果。体现在代码中为PICTURE_ENHANCE_RATIO.

参考

[1]岳亚伟,《数字图像处理与python实现》

腾讯云开发者社区

腾讯云面向开发者汇聚海量精品云计算使用和开发经验，营造开放的云计算技术生态圈。

更多推荐

python函数编写脚本

1.输入一个用户名判断是否存在#!/usr/bin/env python3## user## shanhaiimport subprocessprint ("请输入用户名")x = input(">>:")def user_input(x):while True:if not x:continueif x == 'q':breakret = subprocess.ge

腾讯云开发者社区

直接通过javascript-obfuscator命令行混淆js代码

腾讯云开发者社区

错误：“Cannot load JDBC driver class 'com.mysql.jdbc.Driver”的解决方法

“Cannot load JDBC driver class 'com.mysql.jdbc.Driver ” 表示没有JDBC连接MySql的驱动包，因此需要手动添加驱动包到WEB-INF目录下的lib目录中。从网上下载驱动包mysql-connector-java-5.1.5.zip，解压缩后将mysql-connector-java-5.1.5-bin.jar放到“**/w...