中值滤波(C++实现)_c++中值滤波

本文翻译自librow——Median filter

中值滤波简介

中值滤波是一类非线性的窗口滤波，它可以很容易消除破坏性的噪声而保护信号的边缘。滤波背后的基本思想就是对于任何的信号(图像)元素，找出与它周围最接近的元素。中值滤波的特性很类似于均值滤波但是更擅长处理椒盐噪声和保护边界。另一方面，中值滤波经常被用于斑点噪声的消除，但是也有更有效的技术，如扩散滤波，尽管比较复杂。为了理解如何实现中值滤波，让我们从窗口开始理解。

滤波窗口

让我们想象一下，你可以试着读一下下面模具中的黑色空心字母。

图1 第一个模具

很容易，这个字母的发音是[t]。好的，让我们再来读一下这个字母，但是不同的是这次是另一个模具(stencil)中的字母。

图2 第二个模具

现在这个读音变了，应该是[ð]。下面让我们再来读一次。

图3 第三个模具

这时你的读音应该听起来是这样的[θ]。

你知道这是为什么吗？用数学语言来说，你正在对一个元素(字母t)做一个操作(读)。而这结果(读音）取决于邻近的元素(紧接着t的字母)。

帮助我们选择邻近元素的模具是所谓的窗口。是的，窗口就是模具(stencil or pattern)，通过它，你可以选择邻近的元素(给定元素周围的元素)来帮助你做决定。窗口的另一个名字叫做掩码(mask)，掩码就是模具(stencil)，它可以挡住其他暂时我们不关注的元素。

在这个例子中，我们操作的元素位于窗口的最左边。但是实际上，它通常位于窗口的中心位置。

让我们来看一下几个窗口的例子，下面是一维的窗口：

图4 大小为5的一维窗口

下面是二维的窗口，

图5 大小为3x3的二维窗口

对于三维窗口，你可以想象一下建筑物。现在我们讨论一下建筑物里的房间，这个房间就像一个三维的窗口，他将整个建筑的空间划分为几个子空间。你可以在卷积图像处理发现三维窗口。

图6 大小为3x3x3的三维窗口

理解中值滤波

现在让我们来看一看如何找到与其他元素最相像的元素。基本的思想很简单：排序让后取中值。现在让我们找一下如图7所示的中值。

图7 取中值

对，就是这样，我们刚刚已经用中值滤波过滤了一个一维信号。让我们重新来一遍，一步一步的写下中值滤波的过程。

中值滤波算法：

1. 把窗口放到元素上
2. 取出窗口中的元素
3. 将窗口中的元素排序
4. 选出这些元素中的中值

现在既然我们已经有了算法，那么就可以用代码实现了，让我们一起来编程吧！

一维中值滤波编程

在这一节中，我们研究的是窗口大小为5的一维中值滤波。我们使用大小为N的一维信号作为输入。第一步是放置窗口，这里我们通过改变开始的索引值来放置窗口。

//   Move window through all elements of the signal
for (int i = 2; i < N - 2; ++i)
1
2

注意，我们是从第三个元素开始，倒数第三个元素结束。问题就是我们无法从第一个元素开始，因为在这种情况下，过滤窗口的左半部分是空的。我们将在之后讨论如何解决这个问题。

第二步就是取出窗口中的元素。

//   Pick up window elements
for (int j = 0; j < 5; ++j)
   window[j] = signal[i - 2 + j];
1
2
3

第三步就是把窗口中的元素排序。但是在这里我们将会使用一个代码优化的技巧。由于我们需要的仅仅是中值，所以我们只要给一半的元素排序就可以了。

//   Order elements (only half of them)
for (int j = 0; j < 3; ++j)
{
   //   Find position of minimum element
   int min = j;
   for (int k = j + 1; k < 5; ++k)
      if (window[k] < window[min])
         min = k;
   //   Put found minimum element in its place
   const element temp = window[j];
   window[j] = window[min];
   window[min] = temp;
}
1
2
3
4
5
6
7
8
9
10
11
12
13

最后一步——取中值。

//   Get result - the middle element
result[i - 2] = window[2];
1
2

最后，让我们把整个算法写成函数的形式,

//   1D MEDIAN FILTER implementation
//     signal - input signal
//     result - output signal
//     N      - length of the signal
void _medianfilter(const element* signal, element* result, int N)
{
   //   Move window through all elements of the signal
   for (int i = 2; i < N - 2; ++i)
   {
      //   Pick up window elements
      element window[5];
      for (int j = 0; j < 5; ++j)
         window[j] = signal[i - 2 + j];
      //   Order elements (only half of them)
      for (int j = 0; j < 3; ++j)
      {
         //   Find position of minimum element
         int min = j;
         for (int k = j + 1; k < 5; ++k)
            if (window[k] < window[min])
               min = k;
         //   Put found minimum element in its place
         const element temp = window[j];
         window[j] = window[min];
         window[min] = temp;
      }
      //   Get result - the middle element
      result[i - 2] = window[2];
   }
}
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

类型element可以这样定义：

typedef double element;
1

处理边界问题

对于所有的窗口过滤器，它们有一个共同的问题。那就是处理边界的问题。如果你把窗口放到第一个或最后一个元素上，那么这个窗口的左边或右边就会是空的。为了填补这个间隙，信号需要被扩展一下。有一个对称扩展信号或图像的好办法，就像这样：

图8 信号扩展

所以在将信号交给我们的中值滤波函数处理之前，信号需要被扩展。让我们来写下这个包装器，它可以做好所有的预处理。

//   1D MEDIAN FILTER wrapper
//     signal - input signal
//     result - output signal
//     N      - length of the signal
void medianfilter(element* signal, element* result, int N)
{
   //   Check arguments
   if (!signal || N < 1)
      return;
   //   Treat special case N = 1
   if (N == 1)
   {
      if (result)
         result[0] = signal[0];
      return;
   }
   //   Allocate memory for signal extension
   element* extension = new element[N + 4];
   //   Check memory allocation
   if (!extension)
      return;
   //   Create signal extension
   memcpy(extension + 2, signal, N * sizeof(element));
   for (int i = 0; i < 2; ++i)
   {
      extension[i] = signal[1 - i];
      extension[N + 2 + i] = signal[N - 1 - i];
   }
   //   Call median filter implementation
   _medianfilter(extension, result ? result : signal, N + 4);
   //   Free memory
   delete[] extension;
}
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

正如你所看到的，我们的信号考虑到了很多实际问题。我们首先检查输入参数：信号不能为空而且长度必须是正的。

//   Check arguments
if (!signal || N < 1)
   return;
1
2
3

第二步我们需要考虑当N为1的情况。这是一种特殊情况，因为被扩展的信号长度至少是2。对于长度为1的信号，滤波结果仍是信号本身。还有一点需要注意，就是当输出参数result为空时，也要保证我们的中值滤波适当的处理。

//   Treat special case N = 1
if (N == 1)
{
   if (result)
      result[0] = signal[0];
   return;
}
1
2
3
4
5
6
7

现在让我们为信号扩展分配内存。

//   Allocate memory for signal extension
element* extension = new element[N + 4];
1
2

然后检查内存是否分配成功。

//   Check memory allocation
if (!extension)
   return;
1
2
3

现在我们来建立扩展。

//   Create signal extension
memcpy(extension + 2, signal, N * sizeof(element));
for (int i = 0; i < 2; ++i)
{
   extension[i] = signal[1 - i];
   extension[N + 2 + i] = signal[N - 1 - i];
}
1
2
3
4
5
6
7

最后，扩展就建立好了，下面就可以开始滤波了。

//   Call median filter implementation
_medianfilter(extension, result ? result : signal, N + 4);
1
2

对了，也不要忘记释放内存。

//   Free memory
delete[] extension;
1
2

由于我们使用了标准库中的内存管理函数，所以我们需要包含相应的头文件。

#include <memory.h>
1

二维中值滤波编程

对于二维的中值滤波，我们处理的是二维信号或图像。原理是一样的，只不过用的是二维窗口。窗口的变化影响的仅仅是元素的选取，其余的还是一样的，无非就是先排序再选取中值。让我们一起来用代码实现二维的中值滤波。对于二维的中值滤波，我们设置的窗口的大小为3x3。

//   2D MEDIAN FILTER implementation
//     image  - input image
//     result - output image
//     N      - width of the image
//     M      - height of the image
void _medianfilter(const element* image, element* result, int N, int M)
{
   //   Move window through all elements of the image
   for (int m = 1; m < M - 1; ++m)
      for (int n = 1; n < N - 1; ++n)
      {
         //   Pick up window elements
         int k = 0;
         element window[9];
         for (int j = m - 1; j < m + 2; ++j)
            for (int i = n - 1; i < n + 2; ++i)
               window[k++] = image[j * N + i];
         //   Order elements (only half of them)
         for (int j = 0; j < 5; ++j)
         {
            //   Find position of minimum element
            int min = j;
            for (int l = j + 1; l < 9; ++l)
            if (window[l] < window[min])
               min = l;
            //   Put found minimum element in its place
            const element temp = window[j];
            window[j] = window[min];
            window[min] = temp;
         }
         //   Get result - the middle element
         result[(m - 1) * (N - 2) + n - 1] = window[4];
      }
}
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

二维情况下的边界处理

正如我们在一维的情况下看到的，我们同样需要扩展对输入的图像。扩展的方法是在图像的上方和下方各添加一行，在图像的左方和右方各添加一列。

图9 图像扩展

下面是我们的扩展图像的包装器函数。

//   2D MEDIAN FILTER wrapper
//     image  - input image
//     result - output image
//     N      - width of the image
//     M      - height of the image
void medianfilter(element* image, element* result, int N, int M)
{
   //   Check arguments
   if (!image || N < 1 || M < 1)
      return;
   //   Allocate memory for signal extension
   element* extension = new element[(N + 2) * (M + 2)];
   //   Check memory allocation
   if (!extension)
      return;
   //   Create image extension
   for (int i = 0; i < M; ++i)
   {
      memcpy(extension + (N + 2) * (i + 1) + 1,
         image + N * i,
         N * sizeof(element));
      extension[(N + 2) * (i + 1)] = image[N * i];
      extension[(N + 2) * (i + 2) - 1] = image[N * (i + 1) - 1];
   }
   //   Fill first line of image extension
   memcpy(extension,
      extension + N + 2,
      (N + 2) * sizeof(element));
   //   Fill last line of image extension
   memcpy(extension + (N + 2) * (M + 1),
      extension + (N + 2) * M,
      (N + 2) * sizeof(element));
   //   Call median filter implementation
   _medianfilter(extension, result ? result : image, N + 2, M + 2);
   //   Free memory
   delete[] extension;
}
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

原文地址：http://www.librow.com/articles/article-1

本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。