图像卷积与滤波

作者：zouxy09
来源：****
原文：https://blog.****.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()）。
 要在频域中对一副图像进行滤波，滤波器的大小和图像的大小必须要匹配，这样两者的相乘才容易。因为一般滤波器的大小比图像要小，所以我们需要拓展我们的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，就是蓝色实线。