opencv笔记4：模板运算和常见滤波操作

time:2015年10月04日 星期日 00时00分27秒

# opencv笔记4：模板运算和常见滤波操作

这一篇主要是学习模板运算，了解各种模板运算的运算过程和分类，理论方面主要参考《图像工程——图像处理》（章毓晋）一书第3章，空域增强：模板操作。同时也有个疑问：此书第四章，频域图像增强，讲了低通滤波和高通滤波，然而这些东西和模板运算中的平滑、锐化操作有什么区别？。。。
以下是正文：

模板运算

首先我们把所有图像看作矩阵。
模板一般是nxn（n通常是3、5、7、9等很小的奇数）的矩阵。模板运算基本思路：将原图像中某个像素的值，作为它本身灰度值和其相邻像素灰度值的函数。模板中有一个锚点（anchor point），通常是矩阵中心点，和原图像中待计算点对应；整个模板对应的区域，就是原图像中像素点的相邻区域。模板也称为核（kernel）。

前面的解释翻译成公式就是：

g(x,y)=function(f(x,y), template)

常见的function操作有卷积和排序两种。卷积可以立即为一个map-reduce过程：元素对应相乘（mapper），乘积累加（reducer）。显然，卷积是一个线性操作。
排序操作也不难理解：模板的锚点和待计算点绑定后，邻域内所有点进行排序操作，将排序结果中符合策略规定的作为结果。一般的排序算法是O(n log n)的，不知道是不是因此有人认为模板排序运算不是线性的。其实通常处理的图像像像素值都是unsigned char类型的，是[0,255]之间的非负整数，显然用桶排序是可以O(n)复杂度内完成排序的，依我看也是一种线性运算。
如果模板排序前，需要对应元素和模板元素相乘，然后将乘积排序，那么这时候乘积可能是浮点数，排序就基本上是O(n log n)了，这确实是非线性操作了。不过目前我没有见到类似的操作，也觉得没有什么实际的用处。

滤波

模板运算的效果，可能让图像变好，也可能让图像变坏。我们当然需要好的那种模板运算了：）利用像素本身及其邻域像素的灰度关系进行增强的方法，被称为滤波，滤波使用到的模板就是滤波器。（注意：滤波器是一个模板矩阵，也就是核kernel，而具体的卷积操作还是排序操作，不是滤波器）

滤波和卷积的区别

卷积是滤波的一种实现方式。卷积是一种具体的运算，虽然它其实也是有点一种抽象的表述；而滤波则是比卷积要抽象的描述。

高频分量和低频分量

先看看频率的本意：(狭义概念)频率是单位时间内完成周期性变化的次数。推广开来，（广义概念）频率就是指一定时间内的变化次数。
频率在信号处理领域大量使用。信号处理中的函数自变量是时间；数字图像处理被看作类似信号处理，只不过这里的函数自变量不再是时间，而是换成了图像矩阵的像素灰度值。
原来在信号处理中，从前一秒到后一秒，信号周期性变化的次数，就是频率；相应地，在数字图像处理中，从一个像素点到相邻的一个像素点，灰度值变化的多少，就是频率。
所谓高频分量，就是频率值高，就是像素之间灰度变化大，这通常对应着图像区域边缘等；而低频分量，就是频率值低，就是像素灰度之间灰度变化小，这通常是图像中稳定的区域，是在一个object的内部，同属于一个superpixel...

总之，这样的理解下，高频分量对应图像边缘等像素变化大的像素点；低频分量对应着图像中稳定的区域。

平滑滤波和锐化滤波

平滑滤波能去除高频分量，而锐化滤波能去除低频分量。这么说还是抽象，具体讲是：平滑滤波去处噪声，锐化滤波强化边缘、细节与周围的对比度。
平滑滤波主要包括：线性平滑滤波（方框滤波、均值滤波、高斯滤波等）、非线性平滑滤波（中值滤波、序统计滤波）。opencv中对应boxblur、blur、gaussianblur函数。
锐化滤波主要包括：线性锐化滤波（拉普拉斯算子、高频提升滤波）、非线性锐化滤波（基于梯度的锐化滤波、最大-最小锐化变换等）

注意 实际上对图像进行二维傅立叶变换得到频谱图，就是图像梯度的分布图，当然频谱图上的各点与图像上各点并不存在一一对应的关系，即使在不移频的情况下也是没有。傅立叶频谱图上我们看到的明暗不一的亮点，实际上图像上某一点与邻域点差异的强弱，即梯度的大小，也即该点的频率的大小（可以这么理解，图像中的低频部分指低梯度的点，高频部分相反）。一般来讲，梯度大则该点的亮度强，否则该点亮度弱。这样通过观察傅立叶变换后的频谱图，也叫功率图，我们首先就可以看出，图像的能量分布，如果频谱图中暗的点数更多，那么实际图像是比较柔和的（因为各点与邻域差异都不大，梯度相对较小），反之，如果频谱图中亮的点数多，那么实际图像一定是尖锐的，边界分明且边界两边像素差异较大的。

平滑滤波

线性平滑滤波

opencv现在有3个线性平滑滤波器：方框滤波、均值滤波、高斯滤波

先说均值滤波。均值滤波就是用指定大小的、元素全为1的模糊核，对原图进行卷积操作（其实，就是原图像中当前位置对应的核大小的区域，各个元素相加），然后除以核的元素个数。也可以理解为，其核为：元素全为1、系数为元素个数的矩阵。

然后是方框滤波。是均值滤波的推广：元素全为1、系数为alpha的矩阵。alpha等于核的个数时，就是均值滤波；否则一般取1（哦，为什么要搞这么大？难道不会超出255麻？）

再看高斯滤波。前面两个滤波的核（也叫模板），元素值都只有一种。如果模板的元素不只一种，就是加权线性滤波了。
高斯滤波是加权线性滤波的一种，准确说是：模板元素的分布符合二次高斯分布。高斯分布其实就是正态分布。因为一般认为噪声的分布都符合高斯分布，那么去噪也用符合高斯分布的模板，效果会比较好。
通过为GaussianBlur函数传入sigmaX,sigmaY,size等参数，函数能生成相应的符合高斯分布的模板。然后对原图像和模板进行卷积操作，就得到滤波后的图像。

非线性平滑滤波

opencv现在有2个非线性平滑滤波：中值滤波和双边滤波

中值滤波：模板限定区域内，取像素灰度的中值（我理解为中位数），作为计算结果。中值滤波的效果是，让与周围像素灰度值的差比较大的像素改取与周围像素值接近的值，消除了信号序列（这里是模板框定范围内像素点灰度值）的孤立点。因为不是简单的取均值，产生的模糊更少些，通常能比均值滤波更好地保持图像的细节。

测试发现，在有白色噪声的图像上，中值滤波和均值滤波的对比效果很明显：

#include <iostream>
#include <opencv2/opencv.hpp>

using namespace std;
using namespace cv;

int main(){
    //【1】定义原图像
    //【2】调用模糊函数
    //【3】显示结果

    //【1】定义原图像
    Mat srcImage = imread("/home/chris/workspace/clion/blur_img1.png");
    Mat meanBlurImage, medianBlurImage;

    //【2】调用模糊函数
    blur(srcImage, meanBlurImage, Size(3,3));
    medianBlur(srcImage, medianBlurImage, 3);


    //【3】显示结果
    imshow("原图", srcImage);
    imshow("均值滤波", meanBlurImage);
    imshow("中值滤波", medianBlurImage);
    while(waitKey(1)!='q'){}
    destroyAllWindows();

    return 0;
}

双边滤波：
双边滤波是一种简单的、非迭代的保边平滑过滤器：能够去除图像噪声，同时很好地保持边界。缺点是比其他过滤器慢。
查看opencv官方文档，目前（opencv3.0）中，双边滤波器函数的实现依然有问题，“This filter does not work inplace.”
看起来有点沮丧，不管了，先了解下原理。

双边滤波同时考虑了空间域和值域的差别：空间域给人一种“出身”的感觉，模板框定了你周围的像素点，这些点不管它们灰度值是多少，你总要按相应权重对待它们（按模板中对应元素值来处理），比如均值滤波是“一视同仁”，高斯滤波是“像冲击波一样从自身衰减”。空间域滤波器的效果是，能去除噪声。
而从值域的角度看，给人一种“看后天努力程度”的感觉：对于模板框定的周围像素点，考虑它们的灰度值，而不去官它们当中的“老幼尊卑”。这方面的代表是alpha-截尾均值滤波器。值域滤波器的效果是，能保留边界效果。

双边滤波综合考虑了空间域和值域，其计算公式中的权重系统，是定义域核与值域核的乘积。

锐化滤波

主要包括：线性锐化滤波（拉普拉斯算子、高频提升滤波）、非线性锐化滤波（基于梯度的锐化滤波、最大-最小锐化变换等）
线性锐化滤波可以借助模板卷积实现。对应积分运算的模板卷积可以平滑图像，反过来对应微分运算的模板卷积可以锐化图像。锐化模板系数的取值，应该在中心为正而周围远离中心处为负。

当然，图像锐化还可以用高通滤波法来做，不过不属于模板操作的范围，现在还不懂，以后再说。

非线性锐化滤波

这次先来看看非线性锐化滤波。所谓线性还是非线性，是从最后的结果来看，计算步骤是不是线性的：虽然拉普拉斯算子是二阶差分得到的，但是结果上开它等同于做线性模板卷积运算；而Sobel算子等一阶差分方法，因为要分别考虑x、y两个方向然后再合并，整个步骤没法简化，所以是非线性的。从推导的角度看，要先看一阶差分操作，也就是非线性的几个算子。

梯度锐化

梯度锐化依然是一种模板算法卷积算法，经过一系列推导，并整理出对应的运算模板，就是我们最终需要的。
图像一般是二维矩阵，因此梯度锐化法一般在x、y方向分别计算出梯度幅值Gx、Gy，然后再合并。
用差分来近似微分.比如水平垂直差分:

f'(x) = f(i,j)-f(i+1,j)
f'(y) = f(i,j)-f(i,j+1)

或者交叉差分：

f'(x) = f(i,j)-f(i+1, j+1)
f'(y) = f(i+1,j ) - f(i, j+1) 

x、y两个方向上差分的结果，可以通过距离公式和在一起，用来表示最后的计算结果，比如采用水平垂直差分+曼哈顿距离公式，有：

g(i,j)=|f(i,j)-f(i+1,j)|+|f(i,j)-f(i,j+1)|

而使用交叉差分+欧几里得距离公式，有：

g(i,j)=sqrt( (f(i,j)-f(i+1,j+1))^2 + (f(i+1,j)-f(i,j+1)^2 )

或者使用交叉差分+曼哈顿距离，有：

g(i,j)=|(f(i,j)-f(i+1,j+1)| + |(f(i+1,j)-f(i,j+1)|

上式就是Roberts算子的梯度幅值公式，对应的模板为：

Sx=[1  0]
   [0 -1]

Sy=[0 -1]
   [1  0]
   
Gx=Sx * A   #卷积操作
Gy=Sy * A   #卷积操作

梯度锐化的改进 上述算法计算出来的，是一个考虑了周遭像素的、差分、然后相加的结果。这个结果拿来取代原来的像素值，不一定合适，我们可以对此加以判断，如果它落在某个阈值范围内（比如大于阈值T），才算作有效，否则仍然取原来的f(i,j)灰度值。

梯度锐化的不足 考虑低频区域中的一个点f(i,j)，其g(i,j)通常是0；再噪声点f(i,j)，显然它对应的g(i,j)不为0，而且还比较大。往往噪声点在梯度锐化算法中，被增强的效果比普通点更大。因此使用此算法前尽量去除噪声。

梯队锐化对应的模板？ 显然上面提到的公式，对应的模板是2x2的。这和通常使用的奇阶方阵不一样，是不实用的。无论是x方向还是y方向，2x2的方阵都仅仅是考虑了半邻域，而不是整个邻域。不过这个思路是可取的。Sobel算子就是考虑了整个邻域的一阶算子，当然还可以用拉普拉斯这样的二阶算子。

Prewitt算子

最简单的、考虑了整个邻域的算子，是Prewitt算子，其卷积模板为：

   [-1 0 1]
Sx=[-1 0 1]
   [-1 0 1]
   
   [ 1  1  1]
Sy=[ 0  0  0]
   [-1 -1 -1]

显然，它仅仅在一个方向上考虑了权值，而另一个方向上则没有考虑。Sobel算法则是其的进一步优化。

Soble和Canny

Sobel算子用来计算图像的导数，目的是获得图像边缘。因此，它常常被用于边缘检测。Canny算法是比Sobel更完整的边缘检测算法，包括了预处理（高斯滤波）、后处理（阈值法去除非边缘点）。
为什么用导数检测边缘？考虑一维图像f(x)，其边缘是一个点。作为边缘点，其左右两侧像素灰度的变化滤肯定不一样（否则，就是真的平滑图像了，哪还有什么边界）。如果计算f(x)的导数，f(x)一定是导数中的极值。推广到二维图像这也是成立的。只不过求导公式使用差分公式来近似代替了。Sobel和前面的锐化梯度的区别在于，它在x、y方向分别使用一个模板进行卷积运算，得到Gx、Gy两个分量，然后用G=sqrt(Gx^2 + Gy^2)得到梯度幅值。

Sobel算子，一般取模板:

   [-1 0 1]
Sx=[-2 0 2]
   [-1 0 1]
   
   [ 1  2  1]
Sy=[ 0  0  0]
   [-1 -2 -1]

Canny是指Canny边缘检测算法，其步骤包括：

1. 滤波：读取灰度图像后，使用高斯滤波做平滑处理（去噪）
2. 增强：用一阶偏导的有限差分来计算梯度的幅值和方向（锐化）。这一步可以使用Roberts、Sobel、Prewitt等算子,或者说，使用类似Sobel滤波器的滤波步骤。
3. 检测：使用阈值法将非边缘点去除，获得真正的边缘点
   如此看来，Canny边缘检测是模板操作的综合应用了，既有平滑处理，也有锐化处理。

在Sobel算子中，容易计算：

Gx(x,y)=Sx * A   #卷积操作
Gy(x,y)=Sy * A   #卷积操作
G(x,y)=sqrt(Gx^2+Gy^2)   #梯度幅值
theta=arctan(Gy(x,y)/Gx(x,y))   #梯度幅角

其中，梯度幅角在Canny最后一步的检测判断中使用到。

P.S.:opencv中的Sobel函数，如果指定的模板规格为3x3，按照前面的模板矩阵，会产生明显的不精确的结果。通过调用Scharr算子的模板，效果会好些：

   [-3  0  3]
Sx=[-10 0 10]
   [-3  0  3]
   
   [-3 -10 -3]
Sy=[ 0  0   0]
   [-3 -10 -3]

线性锐化滤波

前面使用的是一阶差分，现在使用二阶差分，同样能得到用来锐化滤波的操作。只不过很巧妙的是，拉普拉斯算子的推导过程到最后发现，能够等价于一个权值模板卷积操作，所以是一阶操作。

拉普拉斯算子

拉普拉斯算子是二阶微分算子，也用于线性锐化滤波:

g(i,j)=f''(x)+f''(y)

仍然使用差分的方式，容易得到：

f''(x) = 2f(i,j)-f(i+1,j)-f(i-1,j)
f''(y) = 2f(i,j)-f(i,j+1)-f(i,j-1)
g(i,j)=4f(i,j)-f(i+1,j)-f(i-1,j)-f(i,j+1)-f(i,j-1)

对应的模板是：

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

这是只考虑4邻域的情况。如果考虑8邻域，对应的模板是：

-1 -1 -1
-1  8 -1
-1 -1 -1

高频提升滤波

用原始图像见去平滑或模糊图像能得到非锐化掩模，将非锐化掩模加到原始图像上能得到锐化图像。更进一步，可以把原始图像放大A倍后再减去平滑图像：

h(x,y) = A*f(x,y)-g(x,y) = (A-1)*f(x,y)+(f(x,y)-g(x,y)) = (A-1)*f(x,y) + mask(x,y)

其实就是稍微复杂点的线性组合了。

自定义滤波

如果想到了什么新的算法，或者纯粹想试一试碰碰运气，可以自定义模板矩阵，扔给opencv的filter2d函数，就可以看到效果了。

ref

RachealZhang:双边滤波器的原理及实现
浅墨:【OpenCV入门教程之九】 非线性滤波专场：中值滤波、双边滤波
MrMystery:图像增强-图像锐化
奋斗斌斌的专栏:Canny边缘检测算法原理及其VC实现详解(一)