专栏地址：『youcans 的图像处理学习课』
文章目录：『youcans 的图像处理学习课 - 总目录』
【youcans 的图像处理学习课】7. 空间域图像滤波
图像滤波是在尽可能保留图像细节特征的条件下对目标图像的噪声进行抑制，是常用的图像预处理操作。
空间域图像增强的方法很多，各有不同的特点和作用。本节介绍空间域滤波的平滑（低通滤波）和锐化（高通滤波）方法。常用的平滑算法有高斯平滑、均值平滑、中值平滑、双边滤波、导向滤波等；常用的锐化算法有钝化掩蔽、拉普拉斯算子、Sobel梯度算子和 Scharr 梯度算子。
本文提供上述各种算法的完整例程和运行结果，并通过一个应用案例示范综合使用多种图像增强方法。
文章目录
1. 图像的相关与卷积运算
滤波通常是指对图像中特定频率的分量进行过滤或抑制。图像滤波是在尽可能保留图像细节特征的条件下对目标图像的噪声进行抑制，是常用的图像预处理操作。
数据采集都会带有一定的噪声，图像的噪声可以理解为灰度值的随机变化。对图像在空间域存在的随机噪声，可以通过平滑技术进行抑制或去除，称为空间域图像滤波。
频率域滤波是通过傅里叶变换方法实现的，而空间域滤波则是通过相关与卷积运算实现。常用的平滑处理算法有基于二维离散卷积的高斯平滑、均值平滑，基于统计方法的中值平滑，保留边缘信息的双边滤波、导向滤波等。
空间滤波器是由邻域和定义的操作构成的，滤波器规定了滤波时采用的邻域形状及该区域内像素值的处理方法。滤波器也被称为 “核”、“模板”、“窗口”、“掩模”、“算子”，一般在信号处理中称为 “滤波器”，在数学领域称为 “核”。线性滤波器就是指基于线性核的滤波，也就是卷积运算。
1.1 相关与卷积运算
滤波器核是指像素周围某一大小的矩形邻域，也称为模板、滑动窗口。
**相关运算（Correlation operation）**是利用模板对图像进行邻域操作：将滤波器模板的中心移动到待处理的像素点，对模板区域的各点加权相乘后求和。
大小为 m*n 的核（模板） w 与图像 f(x,y) 的相关运算 ( w ⋄ f ) ( x , y ) (w \diamond f)(x,y) (w⋄f)(x,y) 的数学描述为：
(
w
⋄
f
)
(
x
,
y
)
=
∑
s
=
−
a
a
∑
t
=
−
b
b
w
(
s
,
t
)
∗
f
(
x
+
s
,
y
+
t
)
(w \diamond f)(x,y) = \sum_{s=-a}^a \sum_{t=-b}^b w(s,t) * f(x+s,y+t)
(w⋄f)(x,y)=s=−a∑at=−b∑bw(s,t)∗f(x+s,y+t)
相关运算的计算步骤如下：
（1）将模板在图像中逐点移动，模板中心移动到被处理的像素点上；
（2）将模板区域中的各点的系数（权值）与图像的像素值相乘，对乘积求和，即加权求和；
（3）将加权求和结果赋值给模板中心的像素。
注意， “相关运算” 中的 “相关” 不是 “有关的”，而是一种特定的数学运算方式。
**卷积运算（Convolution operation）**也是利用模板对图像进行邻域操作，只是把相关运算的模板旋转了 180度。
大小为 m*n 的核（模板） w 与图像 f(x,y) 的卷积运算
(
w
★
f
)
(
x
,
y
)
(w \bigstar f)(x,y)
(w★f)(x,y) 的数学描述为：
(
w
★
f
)
(
x
,
y
)
=
∑
s
=
−
a
a
∑
t
=
−
b
b
w
(
s
,
t
)
∗
f
(
x
−
s
,
y
−
t
)
(w \bigstar f)(x,y) = \sum_{s=-a}^a \sum_{t=-b}^b w(s,t) * f(x-s,y-t)
(w★f)(x,y)=s=−a∑at=−b∑bw(s,t)∗f(x−s,y−t)
卷积运算符合交换律、结合律和分配律，即：
f
★
g
=
g
★
f
f
★
(
g
★
h
)
=
(
f
★
g
)
★
h
f
★
(
g
+
h
)
=
(
f
★
g
)
+
(
f
★
h
)
f \bigstar g = g \bigstar f \\ f \bigstar (g \bigstar h) = (f \bigstar g) \bigstar h \\ f \bigstar (g + h) = (f \bigstar g) + (f \bigstar h)
f★g=g★ff★(g★h)=(f★g)★hf★(g+h)=(f★g)+(f★h)
（本图片来自 “小黑鸭” 《OpenCV学习+常用函数记录②：图像卷积与滤波》，特此致谢。）
1.2 图像的边界扩充
相关和卷积运算都要对图像的边界点要进行特殊处理，就需要将边界进行适当扩充。
函数说明：
OpenCV 中提供了函数 cv.copyMakeBorder 进行边界扩充方式，也可以为图像设置边框。
cv.copyMakeBorder(src, top, bottom, left, right, borderType[, dst[, value]]) → dst
参数说明：
- src：进行边界扩充的图像
- top, bottom, left, right：上侧、下侧、左侧、右侧边界扩充的的宽度（像素数）
- value：当 borderType 为 BORDER_CONSTANT 时，以常量（value）填充扩充的边界，默认值为 (0,0,0)
- borderType 边界扩充的类型
- cv2.BORDER_REPLICATE：复制，复制最边缘像素进行填充（aa | abcdefg | gg），中值滤波采用复制法
- cv2.BORDER_REFLECT：对称法，以图像边缘为轴进行对称填充（cba| abcdefg | gfe）
- cv2.BORDER_REFLECTT_101：倒映法，以图像最边缘像素为轴进行对称填充（dcb| abcdefg | fed），函数 filter2D, blur, GaussianBlur, bilateralFilter 中默认的边界处理方法
- cv2.BORDER_WRAP：用另一侧元素来填充这一侧的扩充边界（efg| abcdefg | ab）
- cv2.BORDER_CONSTANT：以常数（value）作为像素值进行扩充（vv | abcdefg | vv）
例程 1.65：图像的边界扩充
# 1.65 图像的边界扩充
img = cv2.imread("../images/imgRose1.jpg") # 读取彩色图像(BGR)
top = bottom = left = right = 50
imgReplicate = cv2.copyMakeBorder(img, top, bottom, left, right, cv2.BORDER_REPLICATE)
imgReflect = cv2.copyMakeBorder(img, top, bottom, left, right, cv2.BORDER_REFLECT)
imgReflect101 = cv2.copyMakeBorder(img, top, bottom, left, right, cv2.BORDER_REFLECT_101)
imgWrap = cv2.copyMakeBorder(img, top, bottom, left, right, cv2.BORDER_WRAP)
imgConstant = cv2.copyMakeBorder(img, top, bottom, left, right, cv2.BORDER_CONSTANT, value=(200,200,200))
plt.figure(figsize=(9, 6))
plt.subplot(231), plt.axis([-50,562,-50,562]), plt.title('ORIGINAL'), plt.axis('off')
plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
plt.subplot(232), plt.axis('off'), plt.title('REPLICATE')
plt.imshow(cv2.cvtColor(imgReplicate, cv2.COLOR_BGR2RGB))
plt.subplot(233), plt.axis('off'), plt.title('REFLECT')
plt.imshow(cv2.cvtColor(imgReflect, cv2.COLOR_BGR2RGB))
plt.subplot(234), plt.axis('off'), plt.title('REFLECT_101')
plt.imshow(cv2.cvtColor(imgReflect101, cv2.COLOR_BGR2RGB))
plt.subplot(235), plt.axis('off'), plt.title('WRAP')
plt.imshow(cv2.cvtColor(imgWrap, cv2.COLOR_BGR2RGB))
plt.subplot(236), plt.axis('off'), plt.title('CONSTANT')
plt.imshow(cv2.cvtColor(imgConstant, cv2.COLOR_BGR2RGB))
plt.show()
1.3 Scipy 实现二维离散卷积（sp.convolve2d）
Scipy 中提供了函数 sp.convolve2d 实现二维离散卷积的计算。
对于二维离散卷积的运算，Python的科学计算包Scipy提供了函数实现该功能：
convolve2d(in1, in2, mode="full", boundary="fill", fillvalue=0) → dst
参数说明：
- in1：进行卷积运算的图像，二维数组——只能处理单通道图像，如灰度图像
- in2：卷积操作的模板（卷积核），二维数组
- mode：卷积类型，‘full’、‘valid’、‘same’，默认值为 ‘full’
- boundary：边界扩充方式，‘fill’、‘wrap’、‘symm’，默认值为 ‘fill’
- ‘fill’：以常数（fillvalue）作为像素值进行扩充（vv | abcdefg | vv）
- ‘symm’：对称法，以图像边缘为轴进行对称填充（cba| abcdefg | gfe）
- ‘wrap’：用另一侧元素来填充这一侧的扩充边界（efg| abcdefg | ab）
- fillvalue：当 boundary=‘fill’ 时，以以常数（fillvalue）作为像素值进行扩充
例程 1.66：scipy.signal 实现图像的二维卷积
# 1.66 scipy.signal 实现图像的二维卷积
img = cv2.imread("../images/imgLena.tif", flags=0) # # flags=0 读取为灰度图像
kernel = np.array([[-3-3j,0-10j,+3-3j], [-10+0j,0+0j,+10+0j], [-3+3j,0+10j,+3+3j]]) # Gx + j*Gy
# scipy.signal 实现卷积运算
from scipy import signal
convFull = signal.convolve2d(img, kernel, boundary='symm', mode='full') # full 卷积
convValid = signal.convolve2d(img, kernel, boundary='symm', mode='valid') # valid 卷积
convSame = signal.convolve2d(img, kernel, boundary='symm', mode='same') # same 卷积
print(img.shape, convFull.shape, convValid.shape, convSame.shape) # 输出图像大小有区别
plt.figure(figsize=(9, 6))
plt.subplot(131), plt.axis('off'), plt.title('Original'), plt.axis('off')
plt.imshow(img, cmap='gray', vmin=0, vmax=255)
plt.subplot(132), plt.axis('off'), plt.title('Convolve (full)')
plt.imshow(np.absolute(convFull), cmap='gray', vmin=0, vmax=255)
plt.subplot(133), plt.axis('off'), plt.title('Convolve (same)')
plt.imshow(np.absolute(convSame), cmap='gray', vmin=0, vmax=255)
plt.tight_layout()
plt.show()
注意事项：
-
signal.convolve2d 只能对二维矩阵进行卷积操作，因此只能处理灰度图像。如果需要处理彩色图像，可以分别对每一通道进行卷积操作来实现。
-
signal.convolve2d 选择不同卷积类型 ‘full’、‘valid’、‘same’ 时，图像卷积效果的差别并不明显，但图像尺寸大小有区别，这与不同类型时采用不同的边界处理方式有关。
img.shape: (512, 512)
convFull.shape: (514, 514)
convValid.shape: (510, 510)
convSame.shape: (512, 512)
1.4 cv2 实现二维离散卷积（flip 和 filter2D）
使用 OpenCV 中的 cv.flip 和 cv.filter2D 函数也可以实现图像的卷积运算。
函数 cv.flip 实现围绕轴线翻转二维阵列，将图像沿轴线进行轴对称变换，可以将图像沿水平方向、垂直方向、或水平/垂直方向同时进行翻转。例程 1.38 介绍了图像的翻转（镜像）的使用方法。
函数 cv.filter2D 对图像与核（模板）进行相关计算，与函数 cv.flip 共同实现卷积运算。
函数说明：
cv.filter2D(src, ddepth, kernel[, dst[, anchor[, delta[, borderType]]]]) → dst
参数说明：
- src：卷积处理的输入图像，可以是灰度图像，也可以是多通道的彩色图像
- dst：卷积处理的输出图像，大小和类型与 src 相同
- ddepth：目标图像每个通道的深度（数据类型），ddepth=-1 表示与输入图像的数据类型相同
- kernel：卷积操作的模板（卷积核），二维实型数组
- anchor：卷积核的锚点位置，默认值 (-1, -1) 表示以卷积核的中心为锚点
- delta：输出图像的偏移量，可选项，默认值为 0
- borderType：边界扩充的类型
cv.filter2D 可以处理灰度图像，也可以直接处理彩色图像，不需要对每一色彩通道分别操作。
例程 1.67：cv2 实现图像的二维卷积
# 1.67：cv2 实现图像的二维卷积
img = cv2.imread("../images/imgGaia.tif", flags=0) # # flags=0 读取为灰度图像
kernel = np.array([[-1, -1, -1], [-1, 9, -1], [-1, -1, -1]]) # Gx + j*Gy
kFlip = cv2.flip(kernel, -1) # 将卷积核旋转180度
# 使用函数filter2D算出same卷积
imgConv1 = cv2.filter2D(img, -1, kFlip,
anchor=(0,0), borderType=cv2.BORDER_CONSTANT)
imgConv2 = cv2.filter2D(img, -1, kFlip,
anchor=(0,0), borderType=cv2.BORDER_REFLECT)
plt.figure(figsize=(9, 6))
plt.subplot(131), plt.axis('off'), plt.title('Original'), plt.axis('off')
plt.imshow(img, cmap='gray', vmin=0, vmax=255)
plt.subplot(132), plt.axis('off'), plt.title('cv2.filter2D (BORDER_CONSTANT)')
plt.imshow(np.absolute(imgConv1), cmap='gray', vmin=0, vmax=255)
plt.subplot(133), plt.axis('off'), plt.title('cv2.filter2D (BORDER_REFLECT)')
plt.imshow(np.absolute(imgConv2), cmap='gray', vmin=0, vmax=255)
plt.tight_layout()
plt.show()
1.5 可分离卷积核
如果卷积核 w 可以被分解为两个或多个较小尺寸卷积核 w1、w2…，即： w = w 1 ★ w 2 w = w1 \bigstar w2 w=w1★w2，则成为可分离卷积核。
秩为 1 的矩阵可以分解为一个列向量与一个行向量的乘积，因此秩为 1 的卷积核是可分离卷积核。
可分离卷积核 w 与图像 f 的卷积（same 卷积），等于先用 f 与 w1 卷积，再用 w2 对结果进行卷积：
w
★
f
=
(
w
1
★
w
2
)
★
f
=
w
2
★
(
w
1
★
f
)
=
(
w
1
★
f
)
★
w
2
w \bigstar f = (w_1 \bigstar w_2)\bigstar f = w_2 \bigstar (w_1 \bigstar f) = (w_1 \bigstar f)\bigstar w_2
w★f=(w1★w2)★f=w2★(w1★f)=(w1★f)★w2
随着图像尺寸与卷积核尺寸的增大，用分离的卷积核依次对图像进行卷积操作，可以有效地提高运算速度。因此，在二维图像处理中，经常将一个可分离卷积核分解为一维水平核 kernalX 和一维垂直核 kernalY 的乘积。
函数 sepFilter2D 实现可分离核（模板）对图像进行线性滤波。
函数说明：
cv.sepFilter2D( src, ddepth, kernelX, kernelY[, dst[, anchor[, delta[, borderType]]]]) → dst # OpenCV4
该函数先用一维水平核 kernalX 对图像的行进行滤波，再用一维垂直核 kernalY 对图像的列进行滤波。
参数说明：
- src：卷积处理的输入图像，可以是灰度图像，也可以是多通道的彩色图像
- dst：卷积处理的输出图像，大小和类型与 src 相同
- ddepth：目标图像每个通道的深度（数据类型），ddepth=-1 表示与输入图像的数据类型相同
- kernelX：水平卷积核向量，一维实型数组
- kernelY：垂直卷积核向量，一维实型数组
- anchor：卷积核的锚点位置，默认值 (-1, -1) 表示以卷积核的中心为锚点
- delta：输出图像的偏移量，可选项，默认值为 0
- borderType：边界扩充的类型
例程 1.68：可分离核的卷积操作
# 1.68：可分离核的卷积操作
imgList = list(range(0, 36))
imgTest = np.array(imgList).reshape(6, 6)
# 可分离卷积核: kernXY = kernX * kernY
kernX = np.array([[-1, 3, -1]], np.float32) # (1,3)
kernY = np.transpose(kernX) # (3,1)
kernXY = kernX * kernY
print(kernX.shape, kernY.shape, kernXY.shape)
from scipy import signal
# 二维卷积核直接对图像进行卷积操作
imgConv_XY = signal.convolve2d(imgTest, kernXY, mode='same', boundary='fill')
# 可分离卷积核分解为一维水平核 kernalX 和一维垂直核 kernalY 分别进行卷积操作
imgConv_X = signal.convolve2d(imgTest, kernX, mode='same', boundary='fill')
imgConv_X_Y = signal.convolve2d(imgConv_X, kernY, mode='same', boundary='fill')
print("\n比较 imgConv_XY 与 imgConv_X_Y 是否相等：\t", (imgConv_XY == imgConv_X_Y).all())
print("\nimgConv_kernXY:\n", imgConv_XY)
print("\nimgConv_kernX_kernY:\n", imgConv_X_Y)
运行结果如下：
(1, 3) (3, 1) (3, 3)
比较 imgConv_XY 与 imgConv_X_Y 是否相等： True
imgConv_kernXY:
[[-14. -4. -2. 0. 2. 10.]
[ 11. 7. 8. 9. 10. 23.]
[ 23. 13. 14. 15. 16. 35.]
[ 35. 19. 20. 21. 22. 47.]
[ 47. 25. 26. 27. 28. 59.]
[130. 68. 70. 72. 74. 154.]]
imgConv_kernX_kernY:
[[-14. -4. -2. 0. 2. 10.]
[ 11. 7. 8. 9. 10. 23.]
[ 23. 13. 14. 15. 16. 35.]
[ 35. 19. 20. 21. 22. 47.]
[ 47. 25. 26. 27. 28. 59.]
[130. 68. 70. 72. 74. 154.]]
例程 1.69：可分离核的图像卷积
# 1.69：可分离核的图像卷积
img = cv2.imread("../images/imgGaia.tif", flags=1)
# 可分离卷积核: kernXY = kernX * kernY
kernX = np.array([[-1, 3, -1]], np.float32) # (1,3)
kernY = np.transpose(kernX) # (3,1)
kernXY = kernX * kernY
# (1) 可分离卷积核分解为一维水平核 kernalX 和一维垂直核 kernalY 分步进行卷积操作
imgConvY = cv2.filter2D(img, -1, kernY,
anchor=(0,0), borderType=cv2.BORDER_CONSTANT)
imgConv_X_Y = cv2.filter2D(imgConvY, -1, kernX,
anchor=(0,0), borderType=cv2.BORDER_CONSTANT)
# (2) 二维卷积核 kernXY 直接对图像进行卷积操作
imgConv_XY = cv2.filter2D(img, -1, kernXY,
anchor=(0, 0), borderType=cv2.BORDER_CONSTANT)
# (3) 一维水平核 kernalX 和一维垂直核 kernalY 进行可分离卷积核的卷积操作
imgConvSep_XY = cv2.sepFilter2D(img, -1, kernX, kernY,
anchor=(0,0), borderType=cv2.BORDER_CONSTANT)
print("\n比较 imgConv_XY 与 imgConv_X_Y 是否相等：\t", (imgConv_XY == imgConv_X_Y).all())
print("\n比较 imgConvSep_XY 与 imgConv_XY 是否相等：\t", (imgConvSep_XY == imgConv_XY).all())
plt.figure(figsize=(9, 6))
plt.subplot(131), plt.axis('off'), plt.title('cv2.filter2D(kernXY)')
plt.imshow(cv2.cvtColor(imgConv_XY, cv2.COLOR_BGR2RGB))
plt.subplot(132), plt.axis('off'), plt.title('cv2.filter2D (kernX->kernY)')
plt.imshow(cv2.cvtColor(imgConv_X_Y, cv2.COLOR_BGR2RGB))
plt.subplot(133), plt.axis('off'), plt.title('cv2.sepFilter2D(kernX,kernY)')
plt.imshow(cv2.cvtColor(imgConvSep_XY, cv2.COLOR_BGR2RGB))
plt.tight_layout()
plt.show()
运行结果如下：
比较 imgConv_XY 与 imgConv_X_Y 是否相等： False
比较 imgConvSep_XY 与 imgConv_XY 是否相等： True
2. 空间域平滑滤波（低通滤波）
图像滤波是在尽可能保留图像细节特征的条件下对目标图像的噪声进行抑制，是常用的图像预处理操作。
平滑滤波也称为低通滤波，可以抑制图像中的灰度突变，使图像变得模糊，是低频增强的空间域滤波技术。平滑滤波常用于：
- 模糊图像和图像降噪。
- 在图像重取样前平滑图像以减少混淆
- 减少图像中无关的细节
- 平滑因灰度级不足所导致的图像的伪轮廓
线性空间滤波是指图像与滤波器核（卷积核）进行卷积计算。平滑卷积核与图像的卷积类似于积分运算，对图像的邻域进行加权求和，可以实现空间域平滑滤波。
2.1 低通盒式滤波器
盒式核是最简单的可分离低通滤波器核。盒式核的模板区域中各像素点的系数相同，因此也是可分离核。
盒式滤波器结构简单，便于快速实现和实验。但盒式滤波器对透镜模糊特性的近似能力较差，而且往往会沿垂直方向模糊图像。
OpenCV 提供了 cv.blur 函数和 cv.boxFilter 函数实现盒式滤波器核低通滤波。
函数说明：
cv.blur(src, ksize[, dst[, anchor[, borderType]]]) → dst
函数 cv.blur 使用的滤波器核的表达式为：
K
=
1
k
s
i
z
e
.
w
i
d
t
h
∗
k
s
i
z
e
.
h
e
i
g
h
t
[
1
1
⋯
1
1
1
⋯
1
⋮
⋮
⋮
⋮
1
1
⋯
1
]
K= \frac{1}{ksize.width * ksize.height} \begin{bmatrix} 1 & 1 &\cdots &1\\ 1 & 1 &\cdots &1\\ \vdots &\vdots &\vdots &\vdots\\ 1 & 1 &\cdots &1 \end{bmatrix}
K=ksize.width∗ksize.height1
11⋮111⋮1⋯⋯⋮⋯11⋮1
参数说明：
- src：低通滤波输入图像，可以是灰度图像，也可以是多通道的彩色图像
- dst：低通滤波输出图像，大小和类型与 src 相同
- ksize：模糊核的大小，元组 (width, height)，宽度、高度应设为正奇数
- anchor：卷积核的锚点位置，默认值 (-1, -1)，表示以卷积核的中心为锚点
- borderType：边界扩充的类型
函数说明：
cv.boxFilter(src, ddepth, ksize[, dst[, anchor[, normalize[, borderType]]]]) → dst
函数 cv.blur 使用的滤波器核的表达式为：
K
=
α
[
1
1
⋯
1
1
1
⋯
1
⋮
⋮
⋮
⋮
1
1
⋯
1
]
α
=
{
1
k
s
i
z
e
.
w
i
d
t
h
∗
k
s
i
z
e
.
h
e
i
g
h
t
,
if normalize=True
1
,
if normalize=False
K= \alpha \begin{bmatrix} 1 & 1 &\cdots &1 \\1 & 1 &\cdots &1\\ \vdots &\vdots &\vdots &\vdots\\ 1 & 1 &\cdots &1 \end{bmatrix}\\ \alpha = \begin{cases} \frac{1}{ksize.width * ksize.height}&, \text{if normalize=True}\\ 1 &, \text{if normalize=False}\\ \end{cases}
K=α
11⋮111⋮1⋯⋯⋮⋯11⋮1
α={ksize.width∗ksize.height11,if normalize=True,if normalize=False
显然，当 normalize=True 时，函数 cv.blur() 等价于函数 cv.boxFilter(normalize=True) 。
参数说明：
- src：低通滤波输入图像，可以是灰度图像，也可以是多通道的彩色图像
- dst：低通滤波输出图像，大小和类型与 src 相同
- ddepth：输出图像每个通道的深度（数据类型），ddepth=-1 表示与输入图像的数据类型相同
- ksize：模糊核的大小，元组 (width, height)，宽度、高度应设为正奇数
- anchor：卷积核的锚点位置，默认值 (-1, -1)，表示以卷积核的中心为锚点
- normalize：归一化选项，默认值 normalize=True 时进行归一化，否则不作归一化处理
- borderType：边界扩充的类型
例程 1.70：图像的低通滤波—盒式滤波器
# 1.70：图像的低通滤波 (盒式滤波器核)
img = cv2.imread("../images/Fig0515a.tif", flags=0) # # flags=0 读取为灰度图像
kSize = (5, 5)
kernel1 = np.ones(kSize, np.float32) / (kSize[0]*kSize[1]) # 生成归一化盒式核
imgConv1 = cv2.filter2D(img, -1, kernel1) # cv2.filter2D 方法
imgConv2 = cv2.blur(img, kSize) # cv2.blur 方法
imgConv3 = cv2.boxFilter(img, -1, kSize) # cv2.boxFilter 方法 (默认normalize=True)
print("比较 cv2.filter2D 与 cv2.blur 方法结果相同吗？\t", (imgConv1 == imgConv2).all())
print("比较 cv2.blur 与 cv2.boxFilter 方法结果相同吗？\t", (imgConv2 == imgConv3).all())
kSize = (11, 11)
imgConv11 = cv2.blur(img, kSize) # cv2.blur 方法
plt.figure(figsize=(9, 6))
plt.subplot(131), plt.axis('off'), plt.title("Original")
plt.imshow(img, cmap='gray', vmin=0, vmax=255)
plt.subplot(132), plt.axis('off'), plt.title("cv2.blur (kSize=[5,5])")
plt.imshow(imgConv2, cmap='gray', vmin=0, vmax=255)
plt.subplot(133), plt.axis('off'), plt.title("cv2.blur (kSize=[11,11])")
plt.imshow(imgConv11, cmap='gray', vmin=0, vmax=255)
plt.tight_layout()
plt.show()
运行结果如下：
比较 cv2.filter2D 与 cv2.blur 方法结果相同吗？ True
比较 cv2.blur 与 cv2.boxFilter 方法结果相同吗？ True
2.2 低通高斯滤波器
实际应用中要求卷积核是各向同性的（圆对称），其响应与方向无关。高斯核是唯一可分离的圆对称核，因此非常适合图像处理，对于去除图像中的随机噪声非常有效。
高斯核的数学表达式为：
w
(
s
,
t
)
=
G
(
s
,
t
)
=
1
2
π
σ
2
e
−
r
2
/
2
σ
2
w(s,t) = G(s,t) = \frac{1}{2\pi\sigma^2} e^{- {r^2}/{2\sigma ^2}}
w(s,t)=G(s,t)=2πσ21e−r2/2σ2
两个一维高斯函数 f 和 g 的乘积和卷积的均值与标准差如下：
m
f
×
g
=
m
f
σ
g
2
+
m
g
σ
f
2
σ
g
2
+
σ
f
2
,
σ
f
×
g
=
σ
f
2
∗
σ
g
2
σ
g
2
+
σ
f
2
m
f
⋆
g
=
m
f
+
m
g
,
σ
f
⋆
g
2
=
σ
f
2
+
σ
g
2
\begin{aligned} m_{f \times g} &= \frac{m_f \sigma _g^2 + m_g \sigma _f^2}{\sigma _g^2 + \sigma _f^2} &,\sigma_{f \times g} &= \frac{\sigma _f^2 * \sigma _g^2}{\sigma _g^2 + \sigma _f^2}\\ m_{f \star g} &= m_f + m_g &,\sigma_{f \star g} ^2 &= \sigma _f^2 + \sigma _g^2 \end{aligned}
mf×gmf⋆g=σg2+σf2mfσg2+mgσf2=mf+mg,σf×g,σf⋆g2=σg2+σf2σf2∗σg2=σf2+σg2
OpenCV 提供了 cv.GaussianBlur 函数实现高斯核低通滤波器，cv.getGaussianKernel 函数可以计算一维高斯滤波器的系数。
函数说明：
cv.GaussianBlur(src, ksize, sigmaX[, dst[, sigmaY[, borderType]]]) → dst
cv.getGaussianKernel(ksize, sigma[, ktype]) → retval
参数说明：
- src：低通滤波输入图像，可以是灰度图像，也可以是多通道的彩色图像
- dst：低通滤波输出图像，大小和类型与 src 相同
- ksize：模糊核的大小，元组 (width, height)，宽度、高度应设为正奇数
- sigmaX：x 轴方向的高斯核标准差
- sigmaY：y 轴方向的高斯核标准差，可选项
- borderType：边界扩充的类型
- sigma：高斯核的标准差
- retval：返回值，高斯滤波器的系数
注意事项：
-
sigmaY 缺省时 sigmaY=sigmaX；sigmaY=sigmaX=0 时，由 ksize 自动计算并设置 sigmaY, sigmaX 的值。
-
如 sigma 为负值，由 ksize 自动计算并设置 sigma 的值：
sigma = 0.3*((ksize-1)/2 - 1) + 0.8。
例程 1.71：图像的低通滤波—高斯滤波器
# 1.71：图像的低通滤波 (高斯滤波器核)
img = cv2.imread("../images/imgLena.tif", flags=1)
kSize = (5, 5)
imgGaussBlur1 = cv2.GaussianBlur(img, (5,5), sigmaX=10)
imgGaussBlur2 = cv2.GaussianBlur(img, (11,11), sigmaX=20)
# 计算高斯核
gaussX = cv2.getGaussianKernel(5, 0)
gaussXY = gaussX * gaussX.transpose(1, 0)
print("gaussX:\n", gaussX)
print("gaussXY:\n", gaussXY)
plt.figure(figsize=(9, 6))
plt.subplot(131), plt.axis('off'), plt.title("Original")
plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
plt.subplot(132), plt.axis('off'), plt.title("ksize=5, sigma=10")
plt.imshow(cv2.cvtColor(imgGaussBlur1, cv2.COLOR_BGR2RGB))
plt.subplot(133), plt.axis('off'), plt.title("ksize=11, sigma=20")
plt.imshow(cv2.cvtColor(imgGaussBlur2, cv2.COLOR_BGR2RGB))
plt.tight_layout()
plt.show()
运行结果如下：
gaussX:
[[0.0625]
[0.25 ]
[0.375 ]
[0.25 ]
[0.0625]]
gaussXY:
[[0.00390625 0.015625 0.0234375 0.015625 0.00390625]
[0.015625 0.0625 0.09375 0.0625 0.015625 ]
[0.0234375 0.09375 0.140625 0.09375 0.0234375 ]
[0.015625 0.0625 0.09375 0.0625 0.015625 ]
[0.00390625 0.015625 0.0234375 0.015625 0.00390625]]
2.3 非线性滤波—中值滤波（Median filter）
中值滤波是一种非线性滤波方法，是基于统计排序方法的滤波器 。中值滤波法将像素点的邻域内的所有像素点灰度值的中值作为该像素点的灰度值。
注意中值不是平均值，而是按大小排序的中间值。由于需要排序操作，中值滤波消耗的运算时间很长。
中值滤波处理后像素点的灰度值，可能保持不变，也可能改变为邻域内其它像素点的灰度值。
中值滤波对于消除图像中的椒盐噪声非常有效。椒盐噪声也称为脉冲噪声，是随机出现的白点或者黑点，通常是由于影像讯号受到干扰而产生，如脉冲干扰、图像扫描。
OpenCV 提供了 cv.medianBlur 函数实现中值滤波算法。
函数说明：
cv.medianBlur(src, ksize[, dst]) → dst
参数说明：
- src：输入图像，可以是灰度图像，也可以是多通道的彩色图像
- dst：输出图像，大小和类型与 src 相同
- ksize：模糊核的线性大小，大于 1 的奇数
例程 1.73：图像的非线性滤波—中值滤波器
# 1.73：图像的非线性滤波 (中值滤波器)
img = cv2.imread("../images/Fig0335a.tif", flags=1)
imgMedianBlur1 = cv2.medianBlur(img, 3)
imgMedianBlur2 = cv2.medianBlur(img, 7)
plt.figure(figsize=(9, 6))
plt.subplot(131), plt.axis('off'), plt.title("Original")
plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
plt.subplot(132), plt.axis('off'), plt.title("cv2.medianBlur(size=3)")
plt.imshow(cv2.cvtColor(imgMedianBlur1, cv2.COLOR_BGR2RGB))
plt.subplot(133), plt.axis('off'), plt.title("cv2.medianBlur(size=7)")
plt.imshow(cv2.cvtColor(imgMedianBlur2, cv2.COLOR_BGR2RGB))
plt.tight_layout()
plt.show()
2.4 非线性滤波—双边滤波（Bilateral filter）
双边滤波是一种非线性滤波方法，是结合图像的空间邻近度和像素值相似度的一种折衷处理，同时考虑空域信息和灰度相似性，在去除噪声的同时有效地保持边缘清晰锐利，对于人像处理具有美颜功能。
边缘的灰度变化较大，高斯滤波会明显地模糊边缘，对于高频细节的保护较弱。双边滤波器在空间中也采用高斯滤波器，但增加了一个反映像素强度差异的高斯方差 σ d \sigma_d σd ，在边缘附近离的较远的像素对边缘上的像素值影响很小，从而保证了边缘附近的像素值，实现边缘保存（edge preserving）。双边滤波器的权值是空间临近度权值和像素值相似度权值的乘积，因此输出像素依赖于当前被卷积像素的邻域，又取决于被卷积像素的灰度值和邻域像素的灰度值的差。
双边滤波器核的数学表达式为：
g
(
i
,
j
)
=
∑
f
(
k
,
l
)
w
∑
w
w
=
w
s
∗
w
r
w
s
=
e
−
[
(
i
−
k
)
2
+
(
j
−
l
)
2
]
/
2
σ
s
2
w
r
=
e
−
∥
(
f
(
i
,
j
)
−
f
(
k
,
l
)
∥
2
/
2
σ
r
2
g(i,j) = \frac{\sum f(k,l) w}{\sum w}\\ w = ws * wr\\ ws = e^{- [{(i-k)^2+(j-l)^2}]/{2\sigma _s^2}}\\ wr = e^{- \lVert {(f(i,j)-f(k,l)} \rVert ^2/{2\sigma _r^2}}
g(i,j)=∑w∑f(k,l)ww=ws∗wrws=e−[(i−k)2+(j−l)2]/2σs2wr=e−∥(f(i,j)−f(k,l)∥2/2σr2
双边滤波器对于低频信息的滤波效果较好，但不能干净地过滤彩色图像里的高频噪声。
OpenCV 提供了 cv. bilateralFilter 函数可以实现图像的双边滤波。
函数说明：
cv.bilateralFilter(src, d, sigmaColor, sigmaSpace[, dst[, borderType]]) → dst
参数说明：
- src：输入图像，可以是灰度图像，也可以是多通道的彩色图像
- dst：输出图像，大小和类型与 src 相同
- d：滤波核的像素邻域直径。如 d<=0 ，则由从 sigmaSpace 计算得到。
- sigmaColor：滤波器核在颜色空间的方差，反映产生颜色影响的颜色强度区间的大小
- sigmaSpace：滤波器核在坐标空间的方差，反映产生颜色影响的影响空间的大小
- borderType：边界扩充的类型
例程 1.74：图像的非线性滤波—双边滤波器
# 1.74：图像的非线性滤波—双边滤波器
img = cv2.imread("../images/imgFabricNoise.png", flags=1)
imgBiFilter = cv2.bilateralFilter(img, d=0, sigmaColor=100, sigmaSpace=10)
imgMeanFilter = cv2.pyrMeanShiftFiltering(img, sp=15, sr=20)
plt.figure(figsize=(9, 6))
plt.subplot(131), plt.axis('off'), plt.title("Original")
plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
plt.subplot(132), plt.axis('off'), plt.title("cv2.bilateralFilter")
plt.imshow(cv2.cvtColor(imgBiFilter, cv2.COLOR_BGR2RGB))
plt.subplot(133), plt.axis('off'), plt.title("cv2.pyrMeanShiftFiltering")
plt.imshow(cv2.cvtColor(imgMeanFilter, cv2.COLOR_BGR2RGB))
plt.tight_layout()
plt.show()
2.5 非线性滤波—联合双边滤波（Joint bilateral filter）
联合双边滤波是在双边滤波基础上对相似性权重模板进行优化，对于纹理图像的处理效果较好。
联合双边滤波与双边滤波的区别在于：双边滤波是根据图像中不同位置的灰度值差异建立相似性权重模板，再与距离权重模板相乘计算卷积核；联合滤波是先对图像进行高斯平滑，然后根据高斯平滑图像的灰度值差异建立相似性模板，再计算卷积核。即联合双边滤波是根据原始图像的高斯平滑作为引导图片来建立相似性模板。进一步地，如果使用其它引导图片来建立相似性模板，还可以实现其它功能。
OpenCV 在 ximgproc 模块提供了 cv.ximgproc.jointBilateralFilter 函数实现联合双边滤波算法。
函数说明：
cv.ximgproc.jointBilateralFilter(joint, src, d, sigmaColor, sigmaSpace[, dst[, borderType]]) → dst
参数说明：
- src：输入图像，可以是灰度图像，也可以是多通道的彩色图像
- joint：联合滤波的导向图像，大小和类型与 src 相同
- dst：输出图像，大小和类型与 src 相同
- d：滤波核的像素邻域直径
- sigmaColor：滤波器核在颜色空间的方差，反映产生颜色影响的颜色强度区间的大小
- sigmaSpace：滤波器核在坐标空间的方差，反映产生颜色影响的影响空间的大小
- borderType：边界扩充的类型
注意事项：
- OpenCV3 中的 ximgproc 模块提供联合双边滤波算法。ximgproc 属于扩展模块，因此需要安装扩展包（opencv-contrib-python）提供支持。
- 除了主模块，还引入了contrib，其中的ximgproc模块包括了联合双边滤波的算法。因此如果需要使用opencv的联合双边滤波，需要安装opencv-contrib-python包。
参考例程：
（略）
2.6 非线性滤波—导向滤波（Guided filter）
导向滤波又称引导滤波，通过一张引导图片反映边缘、物体等信息，对输入图像进行滤波处理，使输出图像的内容由输入图像决定，但纹理与引导图片相似。
导向滤波的原理是局部线性模型，在保持双边滤波的优势（有效保持边缘，非迭代计算）的同时计算速度很快。
而克服双边滤波速度慢的缺点，
导向滤波（导向滤波）不仅能实现双边滤波的边缘平滑，而且在检测到边缘附近有很好的表现，可应用在图像增强、HDR压缩、图像抠图及图像去雾等场景。
在进行保持边缘滤波时，可以采用原始图像自身或其预处理后的图像作为导向图片。
OpenCV 在 ximgproc 模块提供了 cv.ximgproc.guidedFilter 函数实现导向滤波算法。
函数说明：
cv.ximgproc_guidedFilter.filter(guide, src, d[, eps[, dDepth]) → dst
参数说明：
- src：输入图像，可以是灰度图像，也可以是多通道的彩色图像
- guide：导向图像，大小和类型与 src 相同
- dst：输出图像，大小和类型与 src 相同
- d：滤波核的像素邻域直径
- eps：规范化参数， eps 的平方类似于双边滤波中的 sigmaColor
- dDepth：输出图片的数据深度
参考例程：
（略）
3. 空间域锐化滤波（高通滤波）
3.1 图像的梯度算子
图像模糊通过平滑（加权平均）来实现，类似于积分运算。图像锐化则通过微分运算（有限差分）实现，使用一阶微分或二阶微分都可以得到图像灰度的变化值。
图像锐化的目的是增强图像的灰度跳变部分，使模糊的图像变得清晰。图像锐化也称为高通滤波，通过和增强高频，衰减和抑制低频。图像锐化常用于电子印刷、医学成像和工业检测。
- 恒定灰度区域，一阶导数为零，二阶导数为零；
- 灰度台阶或斜坡起点区域，一阶导数非零，，二阶导数非零；
- 灰度斜坡区域，一阶导数非零，二阶导数为零。
一阶导数、二阶导数的有限差分公式为：
∂ f ∂ x = f ( x + 1 ) − f ( x ) ∂ 2 f ∂ x 2 = f ( x + 1 ) − 2 f ( x ) + f ( x − 1 ) \begin{aligned} \dfrac{\partial f}{\partial x} &= f(x+1) - f(x) \\ \dfrac{\partial ^2 f}{\partial x ^2} &= f(x+1) - 2f(x) + f(x-1) \end{aligned} ∂x∂f∂x2∂2f=f(x+1)−f(x)=f(x+1)−2f(x)+f(x−1)
在图像处理中，一阶导数用梯度（二维列向量）表示：
∇
f
=
g
r
a
d
(
f
)
=
[
g
x
g
y
]
=
[
∂
f
/
∂
x
∂
f
/
∂
x
]
\nabla f = grad(f)=\begin{bmatrix}g_x\\g_y\end{bmatrix}=\begin{bmatrix}\partial f /\partial x \\\partial f /\partial x \end{bmatrix}
∇f=grad(f)=[gxgy]=[∂f/∂x∂f/∂x]
梯度指出了像素值的最大变化率的方向，经常用于工业检测中产品缺陷检测和自动检测的预处理。
图像梯度提取方法简单直接，能够有效的描述图像的原始状态，因此发展出多种图像梯度算子：Roberts、Prewitt、Sobel、Laplacian、Scharr。
3.2 钝化掩蔽
简单地，从原始图像中减去一幅平滑处理的钝化图像，也可以实现图像锐化效果，称为钝化掩蔽。
令
f
~
(
x
,
y
)
\tilde{f}(x,y)
f~(x,y) 表示平滑图像，则：
g
m
a
s
k
(
x
,
y
)
=
f
(
x
,
y
)
−
f
~
(
x
,
y
)
g
(
x
,
y
)
=
f
(
x
,
y
)
+
k
∗
g
m
a
s
k
(
x
,
y
)
，
k
>
0
\begin{align} g_{mask} (x,y) &= f(x,y) - \tilde{f}(x,y) \\ g(x,y) &= f(x,y) + k * g_{mask}(x,y)， k>0 \end{align}
gmask(x,y)g(x,y)=f(x,y)−f~(x,y)=f(x,y)+k∗gmask(x,y)，k>0
当 k>1 时，实现高提升滤波；当 k=1 时，实现钝化掩蔽；k<1时，减弱钝化掩蔽。
因此，钝化掩蔽的实现过程是：
（1）对原始图像进行平滑处理，得到平滑图像；
（2）从原始图像中减去平滑图像，产生掩蔽模板；
（3）将原始图像与掩蔽模板加权相加，得到钝化掩蔽。
原图减去模糊图的结果为模板，输出图像等于原图加上加权后的模板，当权重为1得到非锐化掩蔽，当权重大于1成为高提升滤波。
钝化掩蔽没有直接计算和使用梯度算子，但减法运算具有微分运算的特征，因此本质上是梯度算法，可以实现锐化滤波的效果。
例程 1.77：图像锐化: 钝化掩蔽
# 1.77：图像锐化: 钝化掩蔽
img = cv2.imread("../images/Fig0338a.tif", flags=0)
# 对原始图像进行平滑，GaussianBlur(img, size, sigmaX)
imgGauss = cv2.GaussianBlur(img, (5,5), sigmaX=5)
imgGaussNorm = cv2.normalize(imgGauss,dst=None,alpha=0,beta=255,norm_type=cv2.NORM_MINMAX)
# 掩蔽模板：从原始图像中减去平滑图像
imgMask = img - imgGaussNorm
passivation1 = img + 0.6 * imgMask # k<1 减弱钝化掩蔽
imgPas1 = cv2.normalize(passivation1, None, 0, 255, cv2.NORM_MINMAX)
passivation2 = img + imgMask # k=1 钝化掩蔽
imgPas2 = cv2.normalize(passivation2, None, 0, 255, cv2.NORM_MINMAX)
passivation3 = img + 2 * imgMask # k>1 高提升滤波
imgPas3 = cv2.normalize(passivation3, None, 0, 255, cv2.NORM_MINMAX)
plt.figure(figsize=(10, 7))
titleList = ["1. Original", "2. GaussSmooth", "3. MaskTemplate",
"4. Passivation(k=0.5)", "5. Passivation(k=1.0)", "6. Passivation(k=2.0)"]
imageList = [img, imgGauss, imgMask, imgPas1, imgPas2, imgPas3]
for i in range(6):
plt.subplot(2,3,i+1), plt.title(titleList[i]), plt.axis('off')
plt.imshow(imageList[i], 'gray', vmin=0, vmax=255)
plt.tight_layout()
plt.show()
钝化掩蔽的图像锐化效果如下图所示，注意当 k>1 时，实现高提升滤波；当 k=1 时，实现钝化掩蔽；而当 k<1时则会减弱钝化掩蔽。
3.3 拉普拉斯卷积核（Laplacian）
各向同性卷积核的响应与方向无关。最简单的各向同性导数算子（卷积核）是拉普拉斯算子（Laplace）：
∇ 2 f = ∂ 2 f ∂ x 2 + ∂ 2 f ∂ y 2 ∂ 2 f ∂ x 2 = f ( x + 1 , y ) − 2 f ( x , y ) + f ( x − 1 , y ) ∂ 2 f ∂ y 2 = f ( x , y + 1 ) − 2 f ( x , y ) + f ( x , y − 1 ) ∇ 2 f ( x , y ) = f ( x + 1 , y ) + f ( x − 1 , y ) + f ( x , y + 1 ) + f ( x , y − 1 ) − 4 f ( x , y ) \begin{aligned} \nabla ^2 f &= \dfrac{\partial ^2 f}{\partial x ^2} + \dfrac{\partial ^2 f}{\partial y ^2} \\ \dfrac{\partial ^2 f}{\partial x ^2} &= f(x+1,y) - 2f(x,y) + f(x-1,y) \\ \dfrac{\partial ^2 f}{\partial y ^2} &= f(x,y+1) - 2f(x,y) + f(x,y-1) \\ \nabla ^2 f(x,y) &= f(x+1,y) + f(x-1,y) + f(x,y+1) + f(x,y-1) - 4f(x,y) \end{aligned} ∇2f∂x2∂2f∂y2∂2f∇2f(x,y)=∂x2∂2f+∂y2∂2f=f(x+1,y)−2f(x,y)+f(x−1,y)=f(x,y+1)−2f(x,y)+f(x,y−1)=f(x+1,y)+f(x−1,y)+f(x,y+1)+f(x,y−1)−4f(x,y)
由此可以得到拉普拉斯核 K1。类似地，考虑对角项后可以得到拉普拉斯核 K2。
K 1 = [ 0 1 0 1 − 4 1 0 1 0 ] , K 2 = [ 1 1 1 1 − 8 1 1 1 1 ] , K 3 = [ 0 − 1 0 − 1 4 − 1 0 − 1 0 ] , K 4 = [ − 1 − 1 − 1 − 1 8 − 1 − 1 − 1 − 1 ] K1= \begin{bmatrix} 0 & 1 &0\\ 1 & -4 &1\\ 0 & 1 &0\\ \end{bmatrix}, \ K2= \begin{bmatrix} 1 & 1 &1\\ 1 & -8 &1\\ 1 & 1 &1\\ \end{bmatrix}, \ K3= \begin{bmatrix} 0 & -1 &0\\ -1 & 4 &-1\\ 0 & -1 &0\\ \end{bmatrix}, \ K4= \begin{bmatrix} -1 & -1 &-1\\ -1 & 8 &-1\\ -1 & -1 &-1\\ \end{bmatrix} K1= 0101−41010 , K2= 1111−81111 , K3= 0−10−14−10−10 , K4= −1−1−1−18−1−1−1−1
Laplace 是导数算子，会突出图像中的急剧灰度变化，抑制灰度缓慢变化区域，往往会产生暗色背景下的灰色边缘和不连续图像。将拉普拉斯图像与原图叠加，可以得到保留锐化效果的图像。
拉普拉斯卷积核很容易通过卷积操作 cv. filter_2d 实现，OpenCV 也提供了拉普拉斯算子 cv.Laplacian 来实现。
函数说明：
cv.Laplacian(src, ddepth[, dst[, ksize[, scale[, delta[, borderType]]]]]) → dst
参数说明：
-
src：输入图像，可以是灰度图像，也可以是多通道的彩色图像
-
ddepth：输出图片的数据深度：
-
dst：输出图像，大小和类型与 src 相同
-
ksize：计算二阶导数滤波器的孔径大小，必须为正奇数，可选项
-
scale：缩放比例因子，可选项，默认值为 1
-
delta：输出图像的偏移量，可选项，默认值为 0
-
borderType：边界扩充的类型，注意不支持对侧填充（BORDER_WRAP）
例程 1.78：图像锐化：Laplacian 算子
# 1.78：图像锐化：拉普拉斯算子 (Laplacian)
img = cv2.imread("../images/Fig0338a.tif", flags=0) # NASA 月球影像图
# 使用函数 filter2D 实现 Laplace 卷积算子
kernLaplace = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]]) # Laplacian kernel
imgLaplace1 = cv2.filter2D(img, -1, kernLaplace, borderType=cv2.BORDER_REFLECT)
# 使用 cv2.Laplacian 实现 Laplace 卷积算子
imgLaplace2 = cv2.Laplacian(img, -1, ksize=3)
imgRecovery = cv2.add(img, imgLaplace2) # 恢复原图像
# 二值化边缘图再卷积
ret, binary = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY | cv2.THRESH_TRIANGLE)
imgLaplace3 = cv2.Laplacian(binary, cv2.CV_64F)
imgLaplace3 = cv2.convertScaleAbs(imgLaplace3)
plt.figure(figsize=(9, 6))
plt.subplot(131), plt.axis('off'), plt.title("Original")
plt.imshow(img, cmap='gray', vmin=0, vmax=255)
plt.subplot(132), plt.axis('off'), plt.title("cv.Laplacian")
plt.imshow(imgLaplace2, cmap='gray', vmin=0, vmax=255)
plt.subplot(133), plt.axis('off'), plt.title("thresh-Laplacian")
plt.imshow(imgLaplace3, cmap='gray', vmin=0, vmax=255)
plt.tight_layout()
plt.show()
由于拉普拉斯卷积核很敏感，可以先进行阈值化处理，再进行拉普拉斯卷积。例程对比了直接进行拉普拉斯卷积，与阈值化处理后进行拉普拉斯卷积，结果如下图所示。
3.4 Sobel 梯度算子
Sobel 算子是一种离散的微分算子，是高斯平滑和微分求导的联合运算，抗噪声能力强。
Sobel 梯度算子利用局部差分寻找边缘，计算得到梯度的近似值。先计算水平、垂直方向的梯度 G x = k x ∗ s r c G_x = k_x * src Gx=kx∗src， G y = k y ∗ s r c G_y = k_y * src Gy=ky∗src，再求总梯度 $ G = \sqrt{G_x2+G_y2}$ 。编程实现时，可以用绝对值近似平方根： G = ∣ G x ∣ + ∣ G y ∣ G = |G_x| + |G_y| G=∣Gx∣+∣Gy∣。
Sobel 梯度算子是由平滑算子和差分算子卷积得到，Sobel 梯度算子的卷积核为：
K
x
=
[
−
1
0
1
−
2
0
2
−
1
0
1
]
,
K
y
=
[
−
1
−
2
−
1
0
0
0
1
2
1
]
K_x = \begin{bmatrix} -1 & 0 &1\\ -2 & 0 &2\\ -1 & 0 &1\\ \end{bmatrix}, \ K_y = \begin{bmatrix} -1 &-2 &-1\\ 0 &0 &0\\ 1 &2 &1\\ \end{bmatrix}
Kx=
−1−2−1000121
, Ky=
−101−202−101
Sobel 梯度算子很容易通过卷积操作 cv.filter2D 实现，OpenCV 也提供了函数 cv.Sobel 实现 Sobel 梯度算子。
函数说明：
cv.Sobel(src, ddepth, dx, dy[, dst[, ksize[, scale[, delta[, borderType]]]]]) → dst
参数说明：
- src：输入图像，灰度图像，不适用彩色图像
- dst：输出图像，大小和类型与 src 相同
- ddepth：输出图片的数据深度，由输入图像的深度进行选择
- dx：x 轴方向导数的阶数，1 或 2
- dy：y 轴方向导数的阶数，1 或 2
- ksize：Sobel 卷积核的大小，可选的取值为：1/3/5/7，ksize=-1 时使用 Scharr 算子运算
- scale：缩放比例因子，可选项，默认值为 1
- delta：输出图像的偏移量，可选项，默认值为 0
- borderType：边界扩充的类型，注意不支持对侧填充（BORDER_WRAP）
注意事项：
- ddepth 的设置比较复杂，而且容易出错。限于篇幅本文不做展开。
此外，为了处理微分运算导致的数据异常（超出 [0，255]），OpenCV 提供了 cv.convertScaleAbs 进行饱和运算（saturate）： d s t = s a t u r a t e ( s r c ∗ α + b e t a ) dst = saturate(src * \alpha + beta) dst=saturate(src∗α+beta)。
函数说明：
cv.convertScaleAbs(src[, alpha[, beta]]) → dst
参数说明：
- src：输入图像，可以是灰度图像，也可以是多通道的彩色图像
- dst：输出图像，大小和类型与 src 相同
- alpha：调节系数，可选项，默认值为 1
- beta：亮度调节，可选项，默认值为 0
例程 1.79：图像锐化：Sobel 算子
# 1.79：图像锐化：Sobel 算子
# img = cv2.imread("../images/Fig0338a.tif", flags=0) # NASA 月球影像图
img = cv2.imread("../images/imgGaia.tif", flags=0)
# 使用函数 filter2D 实现 Sobel 算子
kernSobelX = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]]) # SobelX kernel
kernSobelY = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]]) # SobelY kernel
imgSobelX = cv2.filter2D(img, -1, kernSobelX, borderType=cv2.BORDER_REFLECT)
imgSobelY = cv2.filter2D(img, -1, kernSobelY, borderType=cv2.BORDER_REFLECT)
# 使用 cv2.Sobel 实现 Sobel 算子
SobelX = cv2.Sobel(img, cv2.CV_16S, 1, 0) # 计算 x 轴方向
SobelY = cv2.Sobel(img, cv2.CV_16S, 0, 1) # 计算 y 轴方向
absX = cv2.convertScaleAbs(SobelX) # 转回 uint8
absY = cv2.convertScaleAbs(SobelY) # 转回 uint8
SobelXY = cv2.addWeighted(absX, 0.5, absY, 0.5, 0) # 用绝对值近似平方根
plt.figure(figsize=(10, 6))
plt.subplot(141), plt.axis('off'), plt.title("Original")
plt.imshow(img, cmap='gray', vmin=0, vmax=255)
plt.subplot(142), plt.axis('off'), plt.title("SobelX")
plt.imshow(SobelX, cmap='gray', vmin=0, vmax=255)
# plt.imshow(imgSobelX, cmap='gray', vmin=0, vmax=255)
plt.subplot(143), plt.axis('off'), plt.title("SobelY")
plt.imshow(SobelY, cmap='gray', vmin=0, vmax=255)
# plt.imshow(imgSobelY, cmap='gray', vmin=0, vmax=255)
plt.subplot(144), plt.axis('off'), plt.title("SobelXY")
plt.imshow(SobelXY, cmap='gray')
plt.tight_layout()
plt.show()
3.5 Scharr 算子
Scharr 算子也称为 Scharr 滤波器，计算 x 轴或 y 轴方向的图像差分。
Scharr 算子是 Soble 算子在 ksize=3 时的优化，与 Soble 的速度相同，且精度更高。Scharr 算子与 Sobel 算子的不同点是在平滑部分，其中心元素占的权重更重，相当于使用较小标准差的高斯函数，也就是更瘦高的模板。
Scharr 算子的卷积核为：
G
x
=
[
−
3
0
3
−
10
0
10
−
3
0
3
]
,
G
y
=
[
−
3
10
−
3
0
0
10
3
10
3
]
G_x = \begin{bmatrix} -3 & 0 &3\\ -10 & 0 &10\\ -3 & 0 &3\\ \end{bmatrix}, \ G_y = \begin{bmatrix} -3 &10 &-3\\ 0 &0 &10\\ 3 &10 &3\\ \end{bmatrix}
Gx=
−3−10−30003103
, Gy=
−30310010−3103
Scharr 算子很容易通过卷积操作 cv.filter2D 实现，OpenCV 也提供了函数 cv.Scharr 实现 Scharr 算子。
函数说明：
cv.Scharr(src, ddepth, dx, dy[, dst[, scale[, delta[, borderType]]]]) → dst
参数说明：
- src：输入图像
- dst：输出图像，大小和类型与 src 相同
- ddepth：输出图片的数据深度，由输入图像的深度进行选择
- dx：x 轴方向导数的阶数
- dy：y 轴方向导数的阶数
- scale：缩放比例因子，可选项，默认值为 1
- delta：输出图像的偏移量，可选项，默认值为 0
- borderType：边界扩充的类型，注意不支持对侧填充（BORDER_WRAP）
例程 1.80：图像锐化：Scharr 算子
# 1.80：图像锐化：Scharr 算子
# img = cv2.imread("../images/Fig0338a.tif", flags=0) # NASA 月球影像图
img = cv2.imread("../images/imgGaia.tif", flags=0)
# 使用函数 filter2D 实现 Scharr 算子
kernScharrX = np.array([[-3, 0, 3], [-10, 0, 10], [-3, 0, 3]]) # ScharrX kernel
kernScharrY = np.array([[-3, 10, -3], [0, 0, 10], [3, 10, 3]]) # ScharrY kernel
# 使用 cv2.Scharr 实现 Scharr 算子
ScharrX = cv2.Scharr(img, cv2.CV_16S, 1, 0) # 计算 x 轴方向
ScharrY = cv2.Scharr(img, cv2.CV_16S, 0, 1) # 计算 y 轴方向
absX = cv2.convertScaleAbs(ScharrX) # 转回 uint8
absY = cv2.convertScaleAbs(ScharrY) # 转回 uint8
ScharrXY = cv2.addWeighted(absX, 0.5, absY, 0.5, 0) # 用绝对值近似平方根
plt.figure(figsize=(10, 6))
plt.subplot(141), plt.axis('off'), plt.title("Original")
plt.imshow(img, cmap='gray', vmin=0, vmax=255)
plt.subplot(142), plt.axis('off'), plt.title("ScharrX")
plt.imshow(ScharrX, cmap='gray', vmin=0, vmax=255)
plt.subplot(143), plt.axis('off'), plt.title("ScharrY")
plt.imshow(ScharrY, cmap='gray', vmin=0, vmax=255)
plt.subplot(144), plt.axis('off'), plt.title("ScharrXY")
plt.imshow(ScharrXY, cmap='gray')
plt.tight_layout()
plt.show()
4. 低通、高通、带阻、带通
图像滤波是在尽可能保留图像细节特征的条件下对目标图像的噪声进行抑制，是常用的图像预处理操作。
空间域和频率域线性滤波器可以分为四类：低通滤波器、高通滤波器、带通滤波器和带阻滤波器，后三类滤波器都可以由低通滤波器构造：
- 低通滤波器： l p ( x , y ) lp(x,y) lp(x,y)
- 高通滤波器： h p ( x , y ) = δ ( x , y ) − l p ( x , y ) hp(x,y) = \delta(x,y) - lp(x,y) hp(x,y)=δ(x,y)−lp(x,y)
- 带通滤波器： b r ( x , y ) = l p 1 ( x , y ) + h p 2 ( x , y ) br(x,y) = lp_1(x,y) + hp_2(x,y) br(x,y)=lp1(x,y)+hp2(x,y)
- 带阻滤波器： b p ( x , y ) = δ ( x , y ) − b r ( x , y ) bp(x,y) = \delta(x,y) - br(x,y) bp(x,y)=δ(x,y)−br(x,y)
这些传递函数都可以由一个低通滤波器传递函数通过线性变换得到。
本案例中建立一维低通滤波器传递函数，进而可以生成空间滤波器核。
以一幅由以下公式描述的同心圆反射板为例，测试滤波方法的特性：
z
(
x
,
y
)
=
[
1
+
c
o
s
(
x
2
+
y
2
)
]
/
2
，
x
,
y
∈
[
−
8.2
,
8.2
]
z(x,y) = [1 + cos(x^2+y2)]/2，\quad x,y \in [-8.2,8.2]
z(x,y)=[1+cos(x2+y2)]/2，x,y∈[−8.2,8.2]
首先，设计一个一维空间低通滤波器函数
x
=
s
i
n
(
x
)
/
x
,
x
∈
[
−
6
π
,
6
π
]
x = sin(x) / x,\quad x \in [-6\pi,6\pi]
x=sin(x)/x,x∈[−6π,6π]
例程 1.81：低通/高通，带阻/带通
# 1.81：低通/高通，带阻/带通
# 同心圆反射板
height, width = 597, 597
m = int((height-1) / 2) # 298
n = int((width-1) / 2) # 298
X = np.linspace(-8.2, 8.2, height) # 长度：597
Y = np.linspace(-8.2, 8.2, width)
x, y = np.meshgrid(X, Y)
circle = 0.5 * (1 + np.cos(x**2 + y**2))
for i in range(circle.shape[0]):
for j in range(circle.shape[1]):
if np.sqrt((i-m)**2 + (j-n)**2) > m:
circle[i,j] = 0
# 一维低通滤波器函数
hFilter, wFilter = 128, 128
mFilter = int((hFilter-1) / 2) # 63
nFilter = int((wFilter-1) / 2) # 63
x = np.linspace(-6*np.pi, 6*np.pi, hFilter) # 长度：128
y = np.linspace(-6*np.pi, 6*np.pi, wFilter) # 长度：128
scale = 1 # 滤波器缩放因子
xFilter = np.sin(x*scale) / x # 128
# 二维低通滤波器函数
X, Y = np.meshgrid(x, y) # 生成x、y网格化数据
Z = np.sin(np.sqrt(X**2+Y**2)) / np.sqrt(X**2+Y**2)
plt.figure(figsize=(10, 3))
plt.subplot(131), plt.axis('off'), plt.title("Concentric circles(597*597)"), plt.imshow(circle, 'gray')
plt.subplot(132), plt.title("1D low-pass filter"), plt.plot(x, xFilter), plt.axis([-20,20,-0.25,1.25])
ax = plt.subplot(133, projection='3d')
surf = ax.plot_surface(X, Y, Z, cmap=plt.get_cmap("rainbow"), linewidth=0, antialiased=False)
ax.set_title("2D low-pass filter")
plt.show()
5. 空间域图像增强技术的综合应用
空间域图像增强的方法很多，各有不同的特点和作用。对于一幅具体图像，往往要根据图像的实际情况，综合使用几种图像增强的方法，以便达到较为理想的结果。
本节以人体骨骼扫描图像（来自G.E.MedicalSystem）为例，要求对图像进行锐化以显示更多的骨骼细节。原始图像的灰度级比较狭窄，噪声含量大，简单使用一种图像增强方法难以达到理想的结果，需要综合应用空间域图像增强技术：首先使用拉普拉斯变换突出细节，然后使用梯度算子增强突出的边缘，再使用低通滤波器降低噪声，以此为模板得到需要的锐化图像，最后用伽马校正调整灰度级的动态范围。具体步骤如下：
（1）拉普拉斯变换，突出原始图像的细节；
（2）原始图像叠加拉普拉斯变换图像，恢复背景特征；
（3）Sobel 梯度算子，增强突出的边缘；
（4）用盒式滤波器平滑梯度图像；
（5）拉普拉斯与平滑梯度相乘得到掩蔽模板；
（6）原始图像与掩蔽模板叠加，得到锐化图像；
（7）Gamma 变换，增大灰度级的动态范围。
例程 1.82 空间域图像增强技术的综合应用
# 1.82 空间域图像增强技术的综合应用
# 原始图像，人体骨骼扫描图像
img = cv2.imread("../images/bonescan.tif", flags=0) # 人体骨骼扫描图像
# 图 2：拉普拉斯变换，突出细节
kernLaplace = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]], np.int8) # Laplacian kernel
# kernLaplaceD = np.array([[1, 1, 1], [1, -8, 1], [1, 1, 1]], np.int8) # Diagonal Laplacian kernel
Laplacian = cv2.filter2D(img, ddepth=-1, kernel=kernLaplace)
imgLaplacian = np.uint8(cv2.normalize(Laplacian, None, 0, 255, cv2.NORM_MINMAX))
# 图 3：原始图像 + 拉普拉斯变换，恢复背景特征
AddLap = img + imgLaplacian
imgAddLap = np.uint8(cv2.normalize(AddLap, None, 0, 255, cv2.NORM_MINMAX))
# 图 4：Sobel 梯度算子，增强突出的边缘
SobelX = cv2.Sobel(img, cv2.CV_16S, 1, 0) # 计算 x 轴方向
SobelY = cv2.Sobel(img, cv2.CV_16S, 0, 1) # 计算 y 轴方向
absX = cv2.convertScaleAbs(SobelX) # 转回 uint8
absY = cv2.convertScaleAbs(SobelY) # 转回 uint8
SobelXY = cv2.addWeighted(absX, 0.5, absY, 0.5, 0) # 用绝对值近似平方根
imgSobel = np.uint8(cv2.normalize(SobelXY, None, 0, 255, cv2.NORM_MINMAX))
# 图 5：用 (5,5) 盒式滤波器平滑梯度图像
kernelBox = np.ones(5, np.float32) / (5 * 5) # 生成归一化盒式核
SobelBox = cv2.filter2D(img, -1, kernelBox) # cv2.filter2D 方法
imgSobelBox = cv2.normalize(SobelBox, None, 0, 255, cv2.NORM_MINMAX)
# 图 6：图2 与 图5 相乘得到模板 mask，突出了强边缘，相对较低了噪声
mask = imgLaplacian * imgSobelBox
imgMask = np.uint8(cv2.normalize(mask, None, 0, 255, cv2.NORM_MINMAX))
# 图7：原始图像与图 6 相加，得到锐化图像，大部分细节更清晰
passivation = img + imgMask * 0.3
imgPassi = np.uint8(cv2.normalize(passivation, None, 0, 255, cv2.NORM_MINMAX))
# 图8: 幂律变换(Gamma 变换)，增大灰度级的动态范围
epsilon = 1e-5 # 非常小的值以防出现除0的情况
# Gamma = np.zeros_like(imgPassi, dtype=np.float)
Gamma = np.power(imgPassi + epsilon, 0.5)
imgGamma = np.uint8(cv2.normalize(Gamma, None, 0, 255, cv2.NORM_MINMAX))
# 绘图
plt.figure(figsize=(10, 7))
titleList = ["1. Original", "2. Laplacian", "3. Original + Laplacian", "4. Sobel",
"5. Sobel Box5", "6. Sobel mask", "7. Passivation", "8. Gamma correction"]
imageList = [img, imgLaplacian, imgAddLap, imgSobel, imgSobelBox, imgMask, imgPassi, imgGamma]
for i in range(8):
plt.subplot(2,4,i+1), plt.title(titleList[i]), plt.axis('off')
plt.imshow(imageList[i], 'gray', vmin=0, vmax=255)
plt.tight_layout()
plt.show()
版权声明：
本文中部分原始图片来自 Rafael C. Gonzalez “Digital Image Processing, 4th.Ed.”，特此致谢。
youcans@xupt 原创作品，转载必须标注原文链接：(https://blog.csdn.net/youcans/article/details/121422046)
Copyright 2022 youcans, XUPT
欢迎关注 『youcans 的 OpenCV 学习课』 系列，持续更新
【youcans 的 OpenCV 学习课】1. 安装与环境配置
【youcans 的 OpenCV 学习课】2. 图像读取与显示
【youcans 的 OpenCV 学习课】3. 图像的创建与修改
【youcans 的 OpenCV 学习课】4. 图像的叠加与混合
【youcans 的 OpenCV 学习课】5. 图像的几何变换
【youcans 的 OpenCV 学习课】6. 灰度变换与直方图处理
【youcans 的 OpenCV 学习课】7. 空间域图像滤波
【youcans 的 OpenCV 学习课】8. 频率域图像滤波（上）
【youcans 的 OpenCV 学习课】9. 频率域图像滤波（下）
5769
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
