高斯滤波

均值滤波是一种线性滤波器，处理思路也很简单，就是将一个窗口区域中的像素计算平均值，然后将窗口中计算得到的均值设置为锚点上的像素值。
该算法有优点在于效率高，思路简单。同样，缺点也很明显，计算均值会将图像中的边缘信息以及特征信息“模糊”掉，会丢失很多特征。

计算均值滤波时可以采用很多优化手段，例如使用积分图的方法对图像进行预处理，处理过后的图像可以通过O(1)的时间复杂度获取窗口区域中的像素和。如果使用并行以及SSE指令集来进行加速，效果会更快。

下面代码使用简单的卷积方案来实现，既然是计算窗口区域中的像素和，即使用如下卷积核即可。图像的边界部分采用padding操作处理。另外，得到的锚点像素值要进行归一化，即除以窗口尺寸大小。
kernel=⎡⎣⎢111111111⎤⎦⎥(3) kernel=\left[\begin{matrix}1 & 1 & 1 \\1 & 1 & 1 \\1 & 1 & 1\end{matrix}\right] \tag{3}
kernel=
⎣
⎡
       
  
1
1
1
       
  
1
1
1
       
  
1
1
1
       
  
⎦
⎤
       
(3)

void MeanFilter(const Mat &src, Mat &dst, int ksize)//均值滤波
{
        CV_Assert(ksize % 2 == 1);

        int *kernel = new int[ksize*ksize];
        for (int i = 0; i < ksize*ksize; i++)
                kernel[i] = 1;
        Mat tmp;
        int len = ksize / 2;
        tmp.create(Size(src.cols + len, src.rows + len), src.type());//添加边框
        dst.create(Size(src.cols, src.rows), src.type());


        int channel = src.channels();
        uchar *ps = src.data;
        uchar *pt = tmp.data;

        for (int row = 0; row < tmp.rows; row++)//添加边框的过程
        {
                for (int col = 0; col < tmp.cols; col++)
                {
                        for (int c = 0; c < channel; c++)
                        {
                                if (row >= len && row < tmp.rows - len && col >= len && col < tmp.cols - len)
                                        pt[(tmp.cols * row + col)*channel + c] = ps[(src.cols * (row - len) + col - len) * channel + c];
                                else
                                        pt[(tmp.cols * row + col)*channel + c] = 0;
                        }
                }
        }


        uchar *pd = dst.data;
        pt = tmp.data;
        for (int row = len; row < tmp.rows - len; row++)//卷积的过程
        {
                for (int col = len; col < tmp.cols - len; col++)
                {
                        for (int c = 0; c < channel; c++)
                        {
                                short t = 0;
                                for (int x = -len; x <= len; x++)
                                {
                                        for (int y = -len; y <= len; y++)
                                        {
                                                t += kernel[(len + x) * ksize + y + len] * pt[((row + x) * tmp.cols + col + y) * channel + c];
                                        }
                                }

                                pd[(dst.cols * (row - len) + col - len) * channel + c] = saturate_cast<ushort> (t/(ksize*ksize));//防止数据溢出ushort是16为数据
                        }
                }
        }
        delete[] kernel;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
高斯滤波

高斯滤波是一种线性滤波，是常用的一种滤波算法，利用二维高斯函数的分布方式来对图像进行平滑。
高斯滤波的优点可以集中在高斯函数的特点上来看
首先，二维高斯函数是旋转对称的，在各个方向上平滑程度相同，不会改变原图像的边缘走向。
第二，高斯函数是单值函数，高斯卷积核的锚点为极值，在所有方向上单调递减，锚点像素不会受到距离锚点较远的像素影响过大，保证了特征点和边缘的特性。
第三，在频域上，滤波过程中不会被高频信号污染。

滤波过程中使用的高斯函数如下：
G(x,y)=12πσ2e&#8722;x2+y22σ2 G(x,y)=\frac{1}{2{\pi}{\sigma}^2}e^{-\frac{x^2+y^2}{2{\sigma}^2}}
G(x,y)=
2πσ
2

1
       
e
&#8722;
2σ
2

x
2
+y
2

       



在构建卷积核后，要将该卷积核归一化处理，也就是将整个高斯卷积核中的值累加，并将卷积核中的每个值除以累加值。
所以，在编写程序的过程中，前面的12πσ2 \frac{1}{2{\pi}{\sigma}^2}
2πσ
2

1
       
部分由于除法会被约去，直接计算e&#8722;x2+y22σ2 e^{-\frac{x^2+y^2}{2{\sigma}^2}}e
&#8722;
2σ
2

x
2
+y
2

       

即可。
另外，高斯的卷积核并不一定都是对称的，即x方向的方差与y方向的方差可能会不同，这里因具体要求而议，程序中使用的计算公式为:
e&#8722;(i&#8722;x)2+(j&#8722;y)22σxσy e^{-\frac{(i-x)^2+(j-y)^2}{2{\sigma}_x{\sigma}_y}}e
&#8722;
2σ
x
       
σ
y
       

(i&#8722;x)
2
+(j&#8722;y)
2

       


其中，x和y表示卷积核的中心，即锚点，i和j表示卷积核中的各个位置。

下面程序使用上面的高斯函数直接运算，性能并不突出，若想提高程序性能可以使用下面两种方法。

第一种，如果了解卷积神经网络的话其实第一个方法很好想，在谷歌的inception 模型中，构建一个3×3的卷积层使用的是1×3与3×1的两个方向上的卷积层来代替的，不但实现了同样的效果，而且还减少了参数。同样，在计算二维高斯核函数时使用同样的原理，将两个方向的卷积核分开计算，即分别计算x方向与y方向的卷积。计算开销会大大减少。

第二种，如果了解数值分析，高斯函数是可以用其它函数逼近的思路来处理。有人提出了这样的一篇论文解决了这个问题
“Recursive implementation of the Gaussian filter”
即使用递归迭代的方法来逼近高斯函数，效率极高。

以上两种方法在配合sse指令以及并行计算以后，效率会有大幅度提升。

//方差可调节
void MyGaussFilter(const Mat &src, Mat &dst, int ksize, double sigmaX, double sigmaY = 0)
{
        CV_Assert(ksize % 2 == 1);

        if (fabs(sigmaY) < 1e-5)
                sigmaY = sigmaX;

        double *kernel = new double[ksize*ksize];
        int center = ksize / 2;
        double sum = 0;
        for (int i = 0; i < ksize; i++)
        {
                for (int j = 0; j < ksize; j++)
                {
                        kernel[i * ksize + j] = exp(-(i - center)*(i - center) / (2 * sigmaX*sigmaX) - (j - center)*(j - center) / (2 * sigmaY*sigmaY));
                        sum += kernel[i*ksize + j];
                }
        }
        for (int i = 0; i < ksize; i++)
        {
                for (int j = 0; j < ksize; j++)
                {
                        kernel[i*ksize + j] /= sum;
                }
        }
        Mat tmp;
        int len = ksize / 2;
        tmp.create(Size(src.cols + len, src.rows + len), src.type());//添加边框
        dst.create(Size(src.cols, src.rows), src.type());


        int channel = src.channels();
        uchar *ps = src.data;
        uchar *pt = tmp.data;

        for (int row = 0; row < tmp.rows; row++)//添加边框的过程
        {
                for (int col = 0; col < tmp.cols; col++)
                {
                        for (int c = 0; c < channel; c++)
                        {
                                if (row >= len && row < tmp.rows - len && col >= len && col < tmp.cols - len)
                                        pt[(tmp.cols * row + col)*channel + c] = ps[(src.cols * (row - len) + col - len) * channel + c];
                                else
                                        pt[(tmp.cols * row + col)*channel + c] = 0;
                        }
                }
        }

        uchar *pd = dst.data;
        pt = tmp.data;
        for (int row = len; row < tmp.rows - len; row++)//卷积的过程
        {
                for (int col = len; col < tmp.cols - len; col++)
                {
                        for (int c = 0; c < channel; c++)
                        {
                                short t = 0;
                                for (int x = -len; x <= len; x++)
                                {
                                        for (int y = -len; y <= len; y++)
                                        {
                                                t += kernel[(len + x) * ksize + y + len] * pt[((row + x) * tmp.cols + col + y) * channel + c];
                                        }
                                }

                                pd[(dst.cols * (row - len) + col - len) * channel + c] = saturate_cast<ushort> (t);//防止数据溢出ushort是16为数据
                        }
                }
        }
        delete[] kernel;
}