图像卷积与滤波

作者：zouxy09
来源：CSDN
原文：https://blog.csdn.net/zouxy09/article/details/49080029
 线性滤波是图像处理的最基本方法。
 使用一个二维的滤波器矩阵（卷积核）对一个待处理的二维图像进行滤波，方法是对于图像的每一个像素点，计算它及邻域像素和滤波器矩阵的对应元素的乘积，然后求和，作为该像素位置的值。这样就完成了滤波过程。
 对图像和滤波矩阵进行逐个元素相乘再求和的操作就叫卷积或者协相关。卷积和协相关的差别是，卷积需要先对滤波矩阵进行180的翻转，如果矩阵是对称的，二者等价。

卷积计算

 对图像处理而言。存在两大类方法：空域处理和频域处理。空域处理是指直接对原始的像素空间进行计算，频率处理是指先对图像变换到频域，再做滤波等处理。

一、空域处理

1.1 2D卷积

 直接2D卷积就是一开始说的那样，对于图像的每一个像素点，计算它的邻域像素和滤波器矩阵的对应元素的乘积，然后加起来，作为该像素位置的值。

直接的实现也称为暴力实现（brute force），因为它严格按照定义来实现，没有任何优化。当然了，在并行实现里面，它也是比较灵活的。另外，也存在一个优化版本，如果我们的kernel是separable可分的，那么就可以得到一个快5倍左右的卷积方法。

1.2 边界处理

 那卷积核遇到图像边缘怎么办？例如图像顶部的像素，它的上面已经没有像素了，那么它的值如何计算？目前有四种主流的处理方法，我们用一维卷积和均值滤波来说明下。
 我们在1D图像中，用每个像素和它的二邻域的平均值来取代它的值。假设我们有个1D的图像I是这样的：

 对非图像边界的像素的操作比较简单。假设我们对I的第四个像素3做局部平均。也就是我们用2,3和7做平均，来取代这个位置的像素值。也就是，平均会产生一副新的图像J，这个图像在相同位置J (4) = (I(3)+I(4)+I(5))/3 = (2+3+7)/3 = 4。同样，我们可以得到J(3) = (I(2)+I(3)+I(4))/3 =(4+2+3)/3 = 3。需要注意的是，新图像的每个像素都取决于旧的图像，在计算J (4)的时候用J (3)是不对的，而是用I(3)，I(4)和I(5)。所以每个像素都是它和它邻域两个像素的平均。平均是线性的操作，因为每个新的像素都是旧像素的线性组合。

 对卷积，也有必须要考虑的情况是，在图像边界的时候，怎么办？J(1)的值应该是什么？它取决于I(0)，I(1)和I(2)。但是我们没有I(0)呀！图像左边没有值了。有四种方式来处理这个问题：

1. 第一种就是想象I是无限长的图像的一部分，除了我们给定值的部分，其他部分的像素值都是0。在这种情况下，I(0)=0。所以J(1) = (I(0) + I(1) + I(2))/3 = (0 + 5 + 4)/3= 3. 同样，J(10) = (I(9)+I(10)+I(11))/3 = (3+ 6 + 0)/3 = 3.
2. 第二种方法也是想象I是无限图像的一部分。但没有指定的部分是用图像边界的值进行拓展。在我们的例子中，因为图像I最左边的值I(1)=5，所以它左边的所有值，我们都认为是5 。而图像右边的所有的值，我们都认为和右边界的值I(10)一样，都是6。这时候J(1) = (I(0) + I(1) + I(2))/3 = (5 + 5 + 4)/3= 14/3. 而J(10) = (I(9)+I(10)+I(11))/3 = (3 + 6 + 6)/3 = 5。
3. 第三种情况就是认为图像是周期性的。也就是I不断的重复。周期就是I的长度。在我们这里，I(0)和I(10)的值就是一样的，I(11)的值和I(1)的值也是一样的。所以J(1) = (I(0) + I(1) + I(2))/3= (I(10) + I(1)+ I(2))/3 = (6 + 5 + 4)/3 = 5 。
4. 最后一种情况就是不管其他地方了。我们觉得I之外的情况是没有定义的，所以没办法使用这些没有定义的值，所以要使用图像I没有定义的值的像素都没办法计算。在这里，J(1)和J(10)都没办法计算，所以输出J会比原图像I要小。

 这四种方法有各自的优缺点。如果我们想象我们使用的图像只是世界的一个小窗口，然后我们需要使用窗口边界外的值，那么一般来说，外面的值和边界上的值是几乎相似的，所以第二种方法可能更说得过去。

二、频域计算

 这个快速实现得益于卷积定理：时域上的卷积等于频域上的乘积。所以将我们的图像和滤波器通过算法变换到频域后，直接将他们相乘，然后再变换回时域（也就是图像的空域）就可以了。

其中o表示矩阵逐元素相乘。
 那用什么方法将空域的图像和滤波器变换到频域？那就是鼎鼎大名的Fast Fourier Transformation 快速傅里叶变换FFT（在OpenCV里面，使用cv::dft()）。在实际应用中，通常都是用FFT来进行高效的DFT离散傅里叶变换（DFT）计算。
　而对于一幅图像，它是二维的信息，且存在空域中，对它进行DFT变换，可以理解为利用DFT处理二维的信号。在对图片进行了二维DFT后，变换的结果需要使用实数图像加虚数图像或是幅度图像加相位图像的形式，在大多数的实际图像处理过程中，对我们有作用的仅仅是幅度图像，因为它包含了原图像中几乎所有的几何信息成分，在频域中，对于一幅图像，高频的部分代表了图像的细节或是纹理等信息，低频的成分代表了图像的轮廓信息，如果我们对一幅细节众多的图像使用低通滤波器，那么它将失去所有的细节信息，这和信号处理中的滤波器原理是相同的。
 要在频域中对一副图像进行滤波，滤波器的大小和图像的大小必须要匹配，这样两者的相乘才容易。因为一般滤波器的大小比图像要小，所以我们需要拓展我们的kernel，让它和图像的大小一致。

 因为CUDA中的FFT实现是周期的，所以kernel的值也要安排成这样，以支持这种周期性。
 为了保证图像边界的像素也可以得到响应输出，我们也需要拓展我们的输入图像。同时，拓展的方式也要支持周期表达。

　给定N长度的I（image）和K（kernel），为了得到线性卷积，我们需要对I和K进行zero padding（0填充）。为什么要补0，因为DFT假定了输入是无限和周期的，周期是N。

　如上图，对于I和K，如果没有padding的话，隐含着会假定I和K是周期的，以他们的长度N为周期。图中本来N长度的I和K都是黑色虚线的部分，然后如果没有padding，隐含着就会在N之外，加上同样的无数个I，如红色虚线部分，加上了一个周期。对K也是这样。如果是zero padding的话，在黑色虚线的其他地方都全是0了，如图中蓝色部分。将I和K卷积，如果没有padding，如黑色虚线，会有红色那部分的artifact。如果有padding，就是蓝色实线。
　在OpenCV中，实现图像DFT变换的方法比较容易，主要思想是利用几个函数的配合，最主要的函数如下：

1.dft()

dft(InputArray src, OutputArray dst, int flags=0, int nonzeroRows=0);

1. InputArray src: 输入图像，可以是实数或虚数
2. OutputArray dst: 输出图像，其大小和类型取决于第三个参数flags
3. int flags = 0: 转换的标识符，有默认值0.enum { DFT_INVERSE=1, DFT_SCALE=2, DFT_ROWS=4, DFT_COMPLEX_OUTPUT=16, DFT_REAL_OUTPUT=32,DCT_INVERSE = DFT_INVERSE, DCT_ROWS=DFT_ROWS };
   其可取的值如下所示：
   DFT_INVERSE: 用一维或二维逆变换取代默认的正向变换
   DFT_SCALE: 缩放比例标识符，根据数据元素个数平均求出其缩放结果，如有N个元素，则输出结果以1/N缩放输出，常与DFT_INVERSE搭配使用。
   DFT_ROWS: 对输入矩阵的每行进行正向或反向的傅里叶变换；此标识符可在处理多种适量的的时候用于减小资源的开销，这些处理常常是三维或高维变换等复杂操作。
   DFT_COMPLEX_OUTPUT: 对一维或二维的实数数组进行正向变换，这样的结果虽然是复数阵列，但拥有复数的共轭对称性（CCS），可以以一个和原数组尺寸大小相同的实数数组进行填充，这是最快的选择也是函数默认的方法。你可能想要得到一个全尺寸的复数数组（像简单光谱分析等等），通过设置标志位可以使函数生成一个全尺寸的复数输出数组。
   DFT_REAL_OUTPUT: 对一维二维复数数组进行逆向变换，这样的结果通常是一个尺寸相同的复数矩阵，但是如果输入矩阵有复数的共轭对称性（比如是一个带有DFT_COMPLEX_OUTPUT标识符的正变换结果），便会输出实数矩阵。
4. int nonzeroRows = 0: 当这个参数不为0，函数会假设只有输入数组（没有设置DFT_INVERSE）的第一行或第一个输出数组（设置了DFT_INVERSE）包含非零值。这样的话函数就可以对其他的行进行更高效的处理节省一些时间，这项技术尤其是在采用DFT计算矩阵卷积时非常有效。

2.getOptimalDFTSize()

int getOptimalDFTSize(int vecsize);

这个函数返回给定的向量尺寸的傅里叶最优尺寸的大小，这是为了可以使用FFT来加速计算DFT，因此需要扩充图像，而至于扩充多少，就要由这个函数来得到。

1. int vecsize: 输入向量尺寸大小(vector size)
   DFT变换在一个向量尺寸上不是一个单调函数，当计算两个数组卷积或对一个数组进行光学分析，它常常会用0扩充一些数组来得到稍微大点的数组以达到比原来数组计算更快的目的。一个尺寸是2阶指数（2,4,8,16,32…）的数组计算速度最快，一个数组尺寸是2、3、5的倍数（例如：300 = 55322）同样有很高的处理效率。
2. getOptimalDFTSize()函数返回大于或等于vecsize的最小数值N，这样尺寸为N的向量进行DFT变换能得到更高的处理效率。在当前N通过p,q,r等一些整数得出N = $2^p\ast 3^q\ast 5^r$.
   这个函数不能直接用于DCT（离散余弦变换）最优尺寸的估计，可以通过getOptimalDFTSize((vecsize+1)/2)*2得到。

3.copyMakeBorder()

void copyMakeBorder( InputArray src, OutputArray dst,int top, int bottom, int left, int right,int borderType, const Scalar& value=Scalar() );

这个函数的作用的扩充已知图像的边界。

1. InputArray src: 输入图像
2. OutputArray dst: 输出图像，与src图像有相同的类型，其尺寸应为Size(src.cols+left+right, src.rows+top+bottom)
3. int类型的top、bottom、left、right: 在图像的四个方向上扩充像素的值
4. int borderType: 边界类型，由borderInterpolate()来定义，常见的取值为BORDER_CONSTANT
5. const Scalar& value = Scalar(): 如果边界类型为BORDER_CONSTANT则表示为边界值

4.magnitude()

void magnitude(InputArray x, InputArray y, OutputArray magnitude);

1. InputArray x: 浮点型数组的x坐标矢量，也就是实部
2. InputArray y: 浮点型数组的y坐标矢量，必须和x尺寸相同
3. OutputArray magnitude: 与x类型和尺寸相同的输出数组

5.normalize()

void normalize(InputArray src,OutoutArray dst,double alpha =1,double beta = 0,int norm_type=NORM_L2,int dtype=-1,InputArray mask=noArray());

1. InputArray src: 输入图像
2. OutputArray dst: 输出图像，尺寸大小和src相同
3. double alpha = 1: range normalization模式的最小值
4. double beta = 0: range normalization模式的最大值，不用于norm normalization(范数归一化)模式
5. int norm_type = NORM_L2: 归一化的类型，主要有 :
   NORM_INF: 归一化数组的C-范数（绝对值的最大值）
   NORM_L1: 归一化数组的L1-范数（绝对值的和）
   NORM_L2: 归一化数组的L2-范数（欧几里得）
   NORM_MINMAX: 数组的数值被平移或缩放到一个指定的范围，线性归一化，一般较常用。
6. int dtype = -1: 当该参数为负数时，输出数组的类型与输入数组的类型相同，否则输出数组与输入数组只是通道数相同，而depth = CV_MAT_DEPTH(dtype)
7. InputArray mask = noArray(): 操作掩膜版，用于指示函数是否仅仅对指定的元素进行操作。