-- TOC --
图像滤波算法这个词其实很乏，图像缩放算法有时也被成为滤波，还有中文词汇中的各种图像滤镜。本文总结的滤波，主要是各种卷积（convolution）算法，基于opencv。
第一次听到卷积这个词，还是在学习卷积神经网络（CNN）的时候，当时一头雾水。卷积其实就是一个这样的动作：将图像中每一个pixel的值进行修改，修改方式为此pixel以及它周围的pixel的weighted sum。每个pixel周围一圈的pixel都要参与计算，每个位置的权重事先确定（kernel），或者通过学习（CNN）。由于要计算每个pixel周围所有的pixel，中文翻译中就出现了卷这个字。
上图中间那个小正方形，叫做卷积核（convolution kernel，滤波器矩阵），它存放的是每个pixel的权重。注意，中间那个pixel也要参与计算。上图示例中最后得到的数字是-8，如果这是对图像进行卷积，-8这个值可能还要做正规化（normalization）。
Bloody Details of Kernel：
图像滤波的作用：
线性滤波，即每个pixel的值都是周围值的线性计算结果，如果是非线性计算，就是非线性滤波。另一个概念是平移不变形（shift-invariant），即对每一个不同位置的pixel，使用相同的kernel。还有个说法，2D卷积/滤波。
Blurs an image using the normalized box filter. (cv::blur)
这是典型的线性滤波算法，kernel为：
\(Kernel=\cfrac{1}{width\cdot height}\begin{bmatrix}1&1&\cdots&1\\1&1&\cdots&1\\\vdots&\vdots&\ddots&\vdots\\1&1&\cdots&1\end{bmatrix}_{height\times width}\)
width和heigth是kernel的size，注意blur的size可以是偶数。blur就是在kernel覆盖的区域取均值，然后给kernel内的点赋值。
下面是cv::blur接口的python测试代码：
import cv2 as cv
img = cv.imread('11.jpg')
cv.imshow('origin', img)
cv.waitKey(0)
img = cv.blur(img, (5,5))
cv.imshow('blur', img)
cv.waitKey(0)
上面代码中的(5,5)是给kernel设定的size（width,height），越大越模糊。这个size可以不是正方形的，(1,15)可以得到垂直模糊的效果，(15,1)可以得到水平模糊的效果。
border的处理，可以指定borderType这个参数，一般默认的（BORDER_REFLECT_101）就可以了。
下面这张图很nice，说明了各种borderType的含义：
默认的borderType是BORDER_REFLECT_101，它是这样的gfedcb|abcdefgh|gfedcba。
默认的anchor=(-1,-1)，表示kernel的中心点。当kernel为偶数x偶数时，测试发现anchor默认为中心位置的右下那个点。下面这段代码对浮点数进行计算，size取偶数(2,2)，修改anchor的值，观察blur的计算方式：
>>> a
array([[0., 1., 2.],
[3., 4., 5.],
[6., 7., 8.]], dtype=float32)
>>> cv.blur(a, (2,2))
array([[2., 2., 3.],
[2., 2., 3.],
[5., 5., 6.]], dtype=float32)
>>> cv.blur(a, (2,2), anchor=(0,0))
array([[2., 3., 3.],
[5., 6., 6.],
[5., 6., 6.]], dtype=float32)
自己在纸上将border补充完整，然后一个个数字的计算，就能确定每个点是如何计算出来的。使用浮点数是为阻止默认的round操作，以便观察细节，结果上面的计算，恰好都能够除净。
Blurs an image using the box filter. (cv::boxFilter)
boxFilter与blur很相似，只是在接口语法上更灵活。它的kernel定义如下：
\(Kernel=\alpha\cdot\begin{bmatrix}1&1&\cdots&1\\1&1&\cdots&1\\\vdots&\vdots&\ddots&\vdots\\1&1&\cdots&1\end{bmatrix}_{height\times width}\)
其中：
\(\alpha = \begin{cases} \cfrac{1}{width\times height} & \text{when } \text{normalize=true} \\1 & \text{otherwise}\end{cases}\)
当normalize=true的时候，就是单纯的标准的blur。
下面是cv::boxFilter接口的python测试代码：
import cv2 as cv
img = cv.imread('11.jpg')
cv.imshow('origin', img)
cv.waitKey(0)
img = cv.boxFilter(img, -1, (5,5))
cv.imshow('boxFilter', img)
cv.waitKey(0)
Blurs an image using a Gaussian filter. (cv::GaussianBlur) The function convolves the source image with the specified Gaussian kernel. In-place filtering is supported.
数学中的正态分布，也常常被称为高斯分布。这里的高斯模糊，理论上也可以说成是正态模糊。
高斯模糊实质上就是一种均值模糊，只是高斯模糊是按照加权平均的，距离越近的点权重越大，距离越远的点权重越小。通俗的讲，高斯滤波就是对整幅图像进行加权平均的过程，每一个像素点的值，都由其本身和邻域内的其他像素值经过加权平均后得到。
kernel通过指定size和高斯分布的标准差\(\sigma\)得到。
高斯滤波中最重要的参数就是高斯分布的标准差\(\sigma\)。它代表着数据的离散程度，如果\(\sigma\)较小，那么生成的kernel中心系数越大，而周围的系数越小，这样对图像的平滑效果就不是很明显；相反，\(\sigma\)较大时，生成的kernel的各个系数相差就不是很大，比较类似于均值模板，对图像的平滑效果就比较明显。
例如，当两个方向都设定\(\sigma=0\)时，opencv自动生成的3x3的kernel示意如下，此时自动算出\(\sigma=0.8\)：
\(Kernel=\cfrac{1}{16}\cdot\begin{bmatrix}1&2&1\\2&4&2\\1&2&1\end{bmatrix}\)
自动生成kernel的算法细节，请参考cv::getGaussianKernel接口，其中说明了\(\sigma\)通过kernel的size计算得到。我们一般使用默认值即可，除非你知道指定的值的具体含义。
网络上流传甚广的一张二维高斯分布示意图：
下面是cv::GaussianBlur接口的python测试代码：
import cv2 as cv
img = cv.imread('11.jpg')
cv.imshow('origin', img)
cv.waitKey(0)
img = cv.GaussianBlur(img, (5,5), 0)
cv.imshow('GaussianBlur', img)
cv.waitKey(0)
高斯噪声
高斯噪声是指概率分布符合高斯分布的一类噪声。
椒盐噪声是出现在随机位置、噪点深度基本固定的噪声；高斯噪声与其相反，是几乎在每个点上都出现、噪点深度随机的噪声。网上有很多人写教程，给图像增加高斯噪声，基于思路就是准备一个生成符合高斯分布的随机数，然后将这些随机数以此叠加到图像的每个像素点上去。
中值滤波是一种典型的非线性滤波，是基于排序统计理论的一种能够有效抑制噪声的非线性信号处理技术，基本思想是用像素点邻域灰度值的中值来代替该像素点的灰度值（直接用一个中间的值，无其它计算），让周围的像素值接近真实的值从而消除孤立的噪声点。该方法在取出脉冲噪声、椒盐噪声的同时能保留图像的边缘细节。这些优良特性是线性滤波所不具备的。
椒盐噪声（salt & pepper noise）是数字图像的一个常见噪声，所谓椒盐，椒就是黑（0为胡椒噪声），盐就是白（255为盐粒噪声）。椒盐噪声就是在图像上随机出现黑色白色的像素。椒盐噪声是由图像传感器，传输信道，解码处理等产生的黑白相间的亮暗点噪声。
去除脉冲干扰及椒盐噪声最常用的算法是中值滤波。大量的实验研究发现，由摄像机拍摄得到的图像受离散的脉冲、椒盐噪声和零均值的高斯噪声的影响较严重。噪声给图像处理带来很多困难，对图像分割、特征提取、图像识别等具有直接影响。因此，实时采集的图像需进行滤波处理。消除图像中的噪声成份叫做图像的平滑化或滤波操作。滤波的目的有两个：一是抽出对象的特征作为图像识别的特征模式；二是为适应计算机处理的要求，消除图像数字化时所混入的噪声。对滤波处理的要求有两条：一是不能损坏图像轮廓及边缘等重要信息；二是使图像清晰，视觉效果好。
中值滤波对消除椒盐噪声非常有效，能够克服线性滤波器带来的图像细节模糊等弊端，能够有效保护图像边缘信息，是非常经典的平滑噪声处理方法。中值滤波对于细节较多的图像不太适用。
下面是对cv::medianBlur接口的python测试代码：
import cv2 as cv
img = cv.imread('11.jpg')
cv.imshow('origin', img)
cv.waitKey(0)
img = cv.medianBlur(img, 7)
cv.imshow('medianBlur', img)
cv.waitKey(0)
kernel size只需要输入一个尺寸，而且必须是基数，越大越模糊。
均值滤波、方框滤波（boxFilter）、高斯滤波，都是线性滤波方式。由于线性滤波的结果是所有像素值的线性组合，因此含有噪声的像素也会被考虑进去，噪声不会被消除，而是以更柔和的方式存在。这时使用非线性滤波效果可能会更好。
cv::bilateralFilter
高斯滤波只考虑了中心点与周边点的空间距离来计算得到权重。首先，对于图像滤波来说，一个通常的intuition是：（自然）图像在空间中变化缓慢，因此相邻的像素点会更相近，但是这个假设在图像的边缘处变得不成立。如果在边缘处也用这种思路来进行滤波的话，则得到的结果必然会模糊掉边缘，这是不合理的，因此考虑再利用像素点的值的大小进行补充，因为边缘两侧的点的像素值差别很大，因此会使得其加权的时候权重具有很大的差别。可以理解成先根据像素值对要用来进行滤波的邻域做一个分割或分类，再给该点所属的类别相对较高的权重，然后进行邻域加权求和，得到最终结果。双边滤波与高斯滤波相比，对于图像的边缘信息能够更好的保留，其原理为一个与空间距离相关的高斯核函数与一个灰度距离相关的高斯函数相乘。
cv::pyrMeanShiftFiltering
我们可以设计一个低通滤波器，去掉图像中的高频噪声，但是往往也会抑制图像的边缘信号，这就是造成图像模糊的原因。
以均值滤波为例，用均值kernel与图像做卷积。大家都知道，在空间域做卷积，相当于在频域做乘积，而均值kernel在频域是没有高频信号的，只有一个常量的分量，所以均值kernel是对图像局部做低通滤波。
高斯滤波也是一种低通滤波器，因为高斯函数经过傅里叶变换后，在频域的分布依然服从高斯分布。所以它对高频信号有很好的滤除效果。
图像增强需要增强图像的细节，而图像的细节往往就是图像中高频的部分，所以增强图像中的高频信号能够达到图像增强的目的。
图像锐化的目的是使模糊的图像变得更加清晰，其主要方式是增强图像的边缘部分，其实就是增强图像中灰度变化剧烈的部分，所以通过增强图像中的高频信号能够增强图像边缘，从而达到图像锐化的目的。从这里可以看出，可以通过提取图像中的高频信号来得到图像的边缘和纹理信息。
Convolves an image with the kernel. (cv::filter2D) 用一个自己定义的kernel去过滤图像。
下面是实现3x3均值滤波的代码：
import cv2 as cv
import numpy as np
img = cv.imread('11.jpg')
cv.imshow('origin', img)
cv.waitKey(0)
kernel = np.ones((3,3))/9
img = cv.filter2D(img, -1, kernel)
cv.imshow('filter2D-3x3blur', img)
cv.waitKey(0)
下面是实现3x3高斯滤波的代码：
import cv2 as cv
import numpy as np
img = cv.imread('11.jpg')
cv.imshow('origin', img)
cv.waitKey(0)
kg = cv.getGaussianKernel(3,0)
kernel = kg @ kg.T
img = cv.filter2D(img, -1, kernel)
cv.imshow('filter2D-3x3gaussian', img)
cv.waitKey(0)
当给filter2D喂一个偶数x偶数的kernel时，使用默认的anchor，其行为与前面介绍cv.blur接口一样。我理解这样的一致，应该可以推广到整个opencv的接口中去。下面的这段代码，请与cv.blur中最后一段测试代码对照：
>>> a
array([[0., 1., 2.],
[3., 4., 5.],
[6., 7., 8.]], dtype=float32)
>>> k = np.array(((1,1),(1,1)))
>>> k
array([[1, 1],
[1, 1]])
>>> cv.filter2D(a, -1, k)/4
array([[2., 2., 3.],
[2., 2., 3.],
[5., 5., 6.]], dtype=float32)
如果改变默认的anchor值：
>>> cv.filter2D(a, -1, k, anchor=(0,0))/4
array([[2., 3., 3.],
[5., 6., 6.],
[5., 6., 6.]], dtype=float32)
当我们看到一个物体时，首先感受到的就是它的边缘。
图像边缘是图像最基本的特征，所谓边缘(Edge)，是指图像局部特性的不连续性，灰度或结构等信息的突变处称之为边缘。例如，灰度级的突变、颜色的突变、纹理结构的突变等。边缘是一个区域的结束，也是另一个区域的开始，利用该特征可以分割图像。
图像的边缘有方向和幅度两种属性。边缘通常可以通过一阶导数或二阶导数检测得到。一阶导数是以最大值作为对应的边缘的位置，而二阶导数则以过零点作为对应边缘的位置。还有一种就是Canny算子，其是在满足一定约束条件下推导出来的边缘检测最优化算子。
边缘检测的做法，一般都是先将图像转成灰度图，然后在应用某个算子计算，再得到另一张灰度图。边缘检测算法主要是基于图像强度的一阶和二阶导数，但导数通常对噪声很敏感，因此需要采用滤波器来过滤噪声，并调用图像增强或阈值化算法进行处理，最后再进行边缘检测。
对于图像而言，梯度就是两个方向的偏导数组合而成。我们将两个邻近pixel的距离差距定义为1,则：
一阶梯度：
\(\cfrac{\partial f(x,y)}{\partial x} = f(x+1,y) - f(x,y)\)
\(\cfrac{\partial f(x,y)}{\partial y} = f(x,y+1) - f(x,y)\)
二阶梯度：
\(\cfrac{\partial^2 f(x,y)}{\partial x^2} = f(x+1,y)' - f(x,y)'\\ = f(x+2,y) - f(x+1,y) - f(x+1,y) + f(x,y)\\ = f(x+2,y) - 2\cdot f(x+1,y) + f(x,y)\)
将中心点移动到\((x,y)\)，可以得到：
\(\cfrac{\partial^2 f}{\partial x^2} = f(x+1,y) - 2\cdot f(x,y) + f(x-1,y)\)
用同样的方式得到y方向的二阶导数：
\(\cfrac{\partial^2 f}{\partial y^2} = f(x,y+1) - 2\cdot f(x,y) + f(x,y-1)\)
Robert算子是一阶微分算子，而且是第一个边缘检测算子，提出者是Lawrence Roberts in 1963。它是一种利用局部差分算子寻找边缘的算子，它采用对角线方向相邻两象素之差近似梯度幅值检测边缘。检测垂直边缘的效果好于斜向边缘，定位精度高，但对噪声敏感，无法抑制噪声的影响。
Roberts算子是两个2x2的kernel，采用的是对角方向相邻的两个像素之差。kernel分为水平方向和垂直方向，如下所示，从其kernel可以看出，Roberts算子能较好的增强正负45度的图像边缘。
\(\cfrac{\partial{f}}{\partial{x}}=\begin{bmatrix}-1&0\\0&1\end{bmatrix}\)
\(\cfrac{\partial{f}}{\partial{y}}=\begin{bmatrix}0&-1\\1&0\end{bmatrix}\)
像素点\((i,j)\)的\(x\)方向的值为：\((i,j)_x = f(i+1,j+1)-f(i,j)\)
像素点\((i,j)\)的\(y\)方向的值为：\((i,j)_y = f(i,j+1)-f(i+1,j)\)
像素点\((i,j)\)的值为（理论）：\(\sqrt{(i,j)_x^2+(i,j)_y^2}\)
为了简化计算，一般采用（实际）：\(\Big|(i,j)_x\Big|+\Big|(i,j)_y\Big|\)
网络上的很多版本的代码，我认为都是错误的，错在cv.filter2D接口的anchor参数，用默认值就不符合Roberts算子的计算规则。我的测试代码如下：
import cv2 as cv
import numpy as np
img = cv.imread('cat.png')
img = cv.resize(img, None, fx=0.7, fy=0.7, interpolation=cv.INTER_CUBIC)
cv.imshow('origin', img)
gimg = cv.cvtColor(img, cv.COLOR_BGR2GRAY)
cv.imshow('gray', gimg)
kx = np.array(((-1,0),(0,1)))
ky = np.array(((0,-1),(1,0)))
x = cv.filter2D(gimg, -1, kx, anchor=(0,0))
y = cv.filter2D(gimg, -1, ky, anchor=(0,0))
cv.imshow('x', x)
cv.imshow('y', y)
rimg = cv.addWeighted(x, 0.5, y, 0.5, 0)
cv.imshow('Roberts Edge', rimg)
cv.waitKey(0)
cv.destroyWindow('gray')
cv.waitKey(0)
运行效果如下：
x图和y图可以清晰地看出它两的互补关系，要合并在一起形成完整的边缘图像，比如猫嘴里的獠牙。
Prewitt算子是一种图像边缘检测的微分算子，其原理是利用特定区域内像素灰度值产生的差分实现边缘检测。由于Prewitt算子采用3x3的kernel对区域内的像素值进行计算，而Robert算子的kernel为2x2。相比Roberts算子，Prewitt算子对噪声有抑制作用，抑制噪声的原理是通d过像素平均，因此噪声较多的图像处理得比较好。但是像素平均相当于对图像的低通滤波，所以Prewitt算子对边缘的定位却不如Roberts算子。
The x-coordinate is defined here as increasing in the "left"-direction, and the y-coordinate is defined as increasing in the "up"-direction.
水平方向算子：\(G_x=\begin{bmatrix}1&0&-1\\1&0&-1\\1&0&-1\end{bmatrix}=\begin{bmatrix}1\\1\\1\end{bmatrix}\begin{bmatrix}1&0&-1\end{bmatrix}\)
垂直方向算子：\(G_y=\begin{bmatrix}1&1&1\\0&0&0\\-1&-1&-1\end{bmatrix}=\begin{bmatrix}1\\0\\-1\end{bmatrix}\begin{bmatrix}1&1&1\end{bmatrix}\)
上面这段介绍，似乎暗示了一个事实，同一个名称的算子，还有一些其他的形式。比如交换Gx的第1列和第3列，交换Gy的第1行和第3行。
两个方向计算后合并在一起，计算与Roberts算子几乎一样，代码略。
Sobel算子结合了高斯平滑和微分求导，该算子在Prewitt算子的基础上增加了权重的概念，认为相邻点的距离远近对当前像素点的影响是不同的，距离越近的像素点对应当前像素的影响越大，从而实现突出图像边缘轮廓。因为Sobel算子结合了高斯平滑和微分求导（差分），因此结果会具有更多的抗噪性。Sobel算子是一种较为常用的边缘检测方法。
水平方向算子：\(G_x=\begin{bmatrix}1&0&-1\\2&0&-2\\1&0&-1\end{bmatrix}=\begin{bmatrix}1\\2\\1\end{bmatrix}\begin{bmatrix}1&0&-1\end{bmatrix}\)
垂直方向算子：\(G_y=\begin{bmatrix}1&2&1\\0&0&0\\-1&-2&-1\end{bmatrix}=\begin{bmatrix}1\\0\\-1\end{bmatrix}\begin{bmatrix}1&2&1\end{bmatrix}\)
cv::Sobel
scharr与sobel算子思想一样，只是卷积核的系数不同，scharr算子提取边界也更加灵敏，能提取到更细小的边界，但请注意，越是灵敏就越是可能误判。
cv::Scharr
二阶微分算子。
cv::Laplacian
常用算子：
\(H1 = \begin{bmatrix}-1&-1&-1\\-1&c&-1\\-1&-1&-1\end{bmatrix}\)
下面这个用于锐化的算子，也很常见：
\(H2 = \begin{bmatrix}0&-1&0\\-1&5&-1\\0&-1&0\end{bmatrix}\)
测试代码：
import cv2 as cv
import numpy as np
img = cv.imread('cat.png')
img = cv.resize(img, None, fx=0.7, fy=0.7, interpolation=cv.INTER_CUBIC)
cv.imshow('origin', img)
gimg = cv.cvtColor(img, cv.COLOR_BGR2GRAY)
gimg = cv.GaussianBlur(gimg, (3,3), 0)
cv.imshow('gray', gimg)
k1 = np.array(((-1,-1,-1),(-1,8,-1),(-1,-1,-1)))
k2 = np.array(((0,-1,0),(-1,5,-1),(0,-1,0)))
gimg1 = cv.filter2D(gimg, -1, k1)
gimg2 = cv.filter2D(gimg, -1, k2)
cv.imshow('h1', gimg1)
cv.imshow('h2', gimg2)
cv.waitKey(0)
cv.destroyWindow('origin')
cv.waitKey(0)
本文链接：https://cs.pynote.net/ag/image/202210132/
-- EOF --
-- MORE --