中值滤波

中值滤波

Matlab实现图像滤波（叁）：高效的中值滤波

文章目录

• 背景介绍
• 高效的中值滤波算法
• • Step 1
  • Step 2
  • Step 3
  • Step 4
  • Step 5
  • Step 6
  • Step 7
  • Step 8
  • Step 9
• 算法的Matlab实现
• 结果验证

• 背景介绍

  在Matlab实现图像滤波（贰）中介绍到了有关中值滤波的相关知识。

  

  基本思想是：用一个中值滤波器窗口依次遍历整个图像，在每个位置都计算得到一个中值，然后将这些值拼在一起就得到了中值滤波后的图像。中值滤波可以很好地去除椒盐噪声，如：

  

  参考代码：

  % imgOrg = rgb2gray(imread("peppers.png"));
  imgOrg = imread("peppers.png");
  imgOrg = imnoise(imgOrg, 'salt & pepper'); % Adding salt-pepper noise
  subplot(1,3,1);
  imshow(imgOrg);
  title("Originall img");
  
  % median filter
  imgMedian = medianfilter(imgOrg, 3); % Generally, the size of a filter is an odd num
  subplot(1,3,2);
  imshow(imgMedian);
  title("Median filter img");
  
  % Using Matlab API
  subplot(1,3,3);
  [~, ~, c] = size(imgOrg);
  newImg = zeros(size(imgOrg));
  for i = 1:c
      newImg(:,:,i) = medfilt2(imgOrg(:,:,i));
  end
  
  newImg = uint8(newImg);
  imshow(newImg);
  title("Medfilt img by API");
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function imgMedian = medianfilter(img,filterSize)
  % img: image
  % filterSize: if=5, it means 5*5
  
  % Generate new blank image
  [h, w, c] = size(img);
  padSize = (filterSize-1)/2;
  imgMedianT = zeros([h+2*padSize, w+2*padSize, c]);
  imgMedianT(1+padSize:padSize+h, 1+padSize:padSize+w, :) = img;
  imgMedian = zeros([h, w, c]);
  
  for k = 1:c
      for i = 1:h
          for j = 1:w
              block = imgMedianT(i:i-1+filterSize, j:j-1+filterSize, k);
              imgMedian(i, j, k) = median(block, 'all');
          end
      end
  end
  imgMedian = uint8(imgMedian);
  end
  

  Elapsed time is 0.895780 seconds.

  从上面可知，每次移动滤波器窗口就要重新求一次中值，而要求中值需要进行排序，所以开销很大，那么有没有更快地方法呢？请看下文。

  高效的中值滤波算法

  不难发现，每一移动滤波器窗口时，该窗口的变化为：丢掉最左边的一列而新增最右边的一列，对于m*n的中值滤波窗口，有mn-2m个像素没有发生变化，因此不需要重新排序。

  

  由此我们得到新的算法：

  ---

  设滤波器窗口大小为 $x\ast y=3\ast 3$ ，一个图像的像素分布为

  

  则：

  Step 1

  $令t=ceil(\frac{xy}{2})$

  解释： $t=ceil(\frac{3\ast 3}{2})=5$ 为排序后中位数所处的位置。取整可以避免不必要的浮点数运算。

  Step 2

  $将窗口移至一个新行的开始，对其内容排序。建立窗口像素的直方图H，确定其中值m，记下亮度\le m的像素数目n_m$

  解释：当滤波器移至如图示黄色区域窗口时：

  

  $n_m$ 相当于当前中值的位置，例如下图 $3\ast 3$ 滤波器窗口排序好后， $105$ 为中值，此时 $n_m=5$ ，即此时中值的位置为 $5$ 号位置。

  

  Step 3

  $$\begin{gathered} 向右移动滤波器，对于移出的左列的每个亮度为p_g的像素p,做如下操作： \\ H[p_g]=H[p_g]-1，若p_g\le m，则n_m=n_m-1 \end{gathered}$$

  解释：当滤波器移至如图示黄色区域窗口时：

  

  最左侧（如红色部分）一列被移除，最右侧一列被移入。所以执行 $H[0]=H[0]-1,H[84]=H[84]-1,H[105]=H[105]-1$ ，即：

  

  由于 $0\le m,84\le m,105\le m,wherem=105$ ，所以现在的 $n_m=5-1-1-1=2$

  Step 4

  $$\begin{gathered} 对于移入的右列的每个亮度为p_g的像素p,做如下操作： \\ H[p_g]=H[p_g]+1，若p_g\le m，则n_m=n_m+1 \end{gathered}$$

  解释：紧接着 $Step3$ ，执行 $H[55]=H[55]+1,H[120]=H[120]+1,H[255]=H[255]+1$ ，即：

  

  由于 $55\le mwherem=105$ ，所以现在的 $n_m=2+1=3$

  Step 5

  $$\begin{gathered} 若n_m=t，则跳至Step8 \\ 若n_m<t，则跳至Step6 \\ 若n_m>t，则跳至Step7 \end{gathered}$$

  解释：$n_m=t$ ，意味着，移动滤波器后，没有改变中值，所以无需进行额外操作。 $n_m\ne t$ ，意味着，移动滤波器后，改变了中值的位置，所以需要进行额外修正操作。这里， $n_m=3,t=5,n_m<t=$ ，故跳至$Step6$ 。其他情况类似，此处不再举例赘述。

  Step 6

  $$\begin{gathered} 重复执行以下操作： \\ m=m+1 \\ {n}_m=n_m+H[m] \\ 直到n_m\ge t，然后跳至Step8 \end{gathered}$$

  解释： $m从105开始往上增，当m=116时，由于H[116]=1（有一个像素的亮度为116），故n_m加1得到n_m=4。在116120之间，不存在亮度为117119的像素，故H[m]=0，n_m无变化。同理当m=120时，n_m=5，此时n_m\ge t=5，故跳至Step8$

  

  $注意，若有两个像素的亮度都是116，则H[116]=2，n_m应该加2$

  Step 7

  $$\begin{gathered} 重复执行以下操作： \\ {n}_m=n_m-H[m] \\ m=m-1 \\ 直到n_m\le t，然后跳至Step8 \end{gathered}$$

  解释：原理同 $Step6$ ，这里不多赘述。

  Step 8

  $$\begin{gathered} 记录m为当前窗口的中值。然后检查滤波窗口的右列是否为图像的右边界： \\ 若是，跳至Step9；反之，跳至Step3 \end{gathered}$$

  解释：中值滤波窗口需要遍历图像的一整行。

  Step 9

  $$\begin{gathered} 检查滤波窗口的底列是否为图像的下边界： \\ 若是，图像遍历完成，算法结束；反之，跳至Step2 \end{gathered}$$

  解释：中值滤波窗口需要遍历图像的每一行，一行滤波完成后换新行重复操作。

  ---

  2 算法的Matlab实现

  主文件main.m：

  img = imnoise(imread("football.jpg"),'salt & pepper',0.1);
  filterSize = 5; % odd
  tic;
  newImg = effMedian(img, filterSize);
  toc;
  subplot(1,2,1);
  imshow(img);
  title("Original Image");
  subplot(1,2,2);
  imshow(newImg);
  title("Filtered Image");
  

  函数文件：

  function [resultImg] = effMedian(inputImg, fSize)
  %EFFMEDIAN Summary of this function goes here
  %   Detailed explanation goes here
  %   inpuImg: input image
  %   fSize: filter size x*x
  
  fSize = [fSize, fSize];
  [ih, iw, ic] = size(inputImg);     % height, width and channel of inputImg                    
  ch = round((fSize(1)-1)/2);
  cw = round((fSize(2)-1)/2);
  
  % padding inputImg
  tempImg = zeros([ih+ch*2, iw+cw*2, ic]);
  tempImg(1+ch:ih+ch, 1+cw:iw+cw, :) = inputImg;
  resultImg = zeros([ih, iw, ic]);
  
  % step1
  t = ceil(fSize(1)*fSize(2)/2);
  
  for c = 1:ic
      for h = 1+ch:ih+ch  % Step 9
          w = 1+cw;
          H = zeros([1,256]); % histogram
          filterBlock = tempImg(h-ch:h+ch, w-cw:w+cw, c);
          tempBlock = filterBlock(:);
          for tmpi = 1:fSize(1)*fSize(2)
              H(tempBlock(tmpi)+1) = H(tempBlock(tmpi)+1)+1;% +1：亮度0对应位置1 
          end
          m = median(tempBlock);
          n_m = length(find(tempBlock <= m));
          resultImg(h-ch,w-1,c) = m;
          for w = 2+cw:iw+cw  % Step 8
              moveOut = tempImg(h-ch:h+ch,w-1-cw,c);
              moveIn = tempImg(h-ch:h+ch,w+cw,c);
              tempBlock = moveOut(:);
              for tmpi = 1:fSize(1)  % Step 3
                  H(tempBlock(tmpi)+1) = H(tempBlock(tmpi)+1)-1;   
              end
              n_m = n_m - length(find(tempBlock <= m)); % Step 3
              tempBlock = moveIn(:);
              for tmpi = 1:fSize(1)  % Step 4
                  H(tempBlock(tmpi)+1) = H(tempBlock(tmpi)+1)+1;   
              end
              n_m = n_m + length(find(tempBlock <= m));  % Step 4
              jmpFlag = n_m - t;                         % Step 5
              if (jmpFlag < 0)                           % Step 6
                  while(n_m < t)
                      m = m + 1;
                      n_m = n_m + H(m+1);
                  end
              elseif (jmpFlag > 0)                       % Step 7
                  while(n_m > t)
                      n_m = n_m - H(m+1);
                      m = m - 1;
                  end
              end
              resultImg(h-cw,w-cw,c) = m;
          end
      end
  end
  
  resultImg = uint8(resultImg);
  
  end
  

  3 结果验证

  

  Elapsed time is 0.759880 seconds.

  显然，高效的中值滤波效果好的同时，实现了更少的开销。

  4 工程下载地址

  见文内代码