自適應中值濾波的原理

     自適應中值濾波的思想是根據噪聲密度改變濾波視窗的大小，同時對噪聲點和訊號點採取不同的處理方法。對噪聲點進行中值濾波，對訊號點保持其灰度值不變。

       設為fij為點(i,j)的灰度值，Sij為當前工作視窗,fmin,fmax和fmed分別為Sij中的灰度最小值、灰度最大值和灰度中值，令maxize為預設的允許最大視窗。自適應中值濾波的步驟如下：

1）若 fmin< fmed <fmax，則轉至第2步；否則增大視窗的尺寸。若的尺寸小於的尺寸，則重複第1步；否則輸出。

2）若 fmin< fij <fmax，則輸出fij；否則輸出fmed

可以看出，演算法中噪聲的檢測和認定時以fmin和fmax為基準的，如果 fmin< fmed <fmax，表明fmed 不是噪聲，

接著根據fmin< fij <fmax，判斷fij 是否為噪聲，當fmed與fij 都不是脈衝噪聲時，優先輸出fij。

引入自適應中值濾波演算法主要有3 個目的: 

一是去除脈衝噪聲；

二是平滑其他非脈衝噪聲；

三是減少諸如物體邊界細化或粗化等失真。

自適應中值濾波的流程圖如下圖所示。

參考程式碼：

/////////////自適應中值濾波/////////////////////////////////
int adp_media_filter(unsigned char* inbuffer,int width,int height,int maxwinsize,unsigned char* outbuffer)
{
	int pos = (maxwinsize - 1) / 2;
	memcpy(outbuffer,inbuffer,width*height);       
	for (int m = pos; m < height - pos; m++)//當前中心位置（m,n）
	{
		for (int n = pos; n < width - pos; n++)
		{
			int curwinsize = 3;     //設定初始濾波視窗大小
			while (curwinsize <= maxwinsize)
			{
				int curpos = (curwinsize - 1) / 2;
				int winpos = 0;
				int lens = curwinsize*curwinsize;
				int* windows = new int[lens];
				for (int i = -curpos; i < curpos + 1; i++)
					for (int j = -curpos; j < curpos + 1; j++)
						windows[winpos++] = inbuffer[(m + i)*width + n + j];
				
				sort(windows, lens);
				int fmin = windows[0];
				int fmax = windows[lens - 1];
				int fmed = windows[(lens - 1) / 2];
				int A1 = fmed - fmin;
				int A2 = fmed - fmax;
				if (A1 > 0 && A2 < 0)//第一層噪聲檢測，fmed是不是噪聲
				{
					//滿足fmin< fmed < fmax，表明fmed不是噪聲
					int B1 = inbuffer[m*width + n] - fmin;//當前視窗中心值inbuffer[m*width + n]
					int B2 = inbuffer[m*width + n] - fmax;
					//滿足fmin< fmn < fmax，表明fmn不是噪聲
					if (B1 > 0 && B2 < 0)//第二層檢測，fmn是不是噪聲
						outbuffer[m*width + n] = inbuffer[m*width + n];//fmn和fmed都不是噪聲優先輸出當前視窗中心值
					else
						outbuffer[m*width + n] = fmed;//fmn是噪聲，輸出中間值重新估計該點
					delete[] windows;
					windows = NULL;
					break;
				}
				curwinsize += 2;
				delete[] windows;
				windows = NULL;
			}
		}
	}
	//對邊界進行處理，與中值濾波一樣
	for (int k = 0;k < pos;k++)
		for (int l =pos;l < width-pos;l++)
			outbuffer[k*width+l] = outbuffer[pos*width+l];
	for(int a=height-pos;a < height;a++)
		for(int b=pos;b < width-pos;b++)
			outbuffer[a*width+b] = outbuffer[(height-pos-1)*width+b];
	for(int c = 0;c < pos;c++)
		for(int d=0;d < height;d++)
			outbuffer[d*width+c] = outbuffer[d*width+pos];
	for (int e = width-pos;e < width;e++)
		for(int f = 0;f < height;f++)
			outbuffer[f*width+e] = outbuffer[f*width+width-pos-1];

	return 0;	
}
 

實驗結果：

2，對乘性噪聲影象進行處理的結果（左為噪聲圖）

注意該噪聲為matlab中的乘性噪聲，可以發現對乘性噪聲的處理效果相當差。

3，對椒鹽噪聲影象進行處理的結果（左為噪聲圖）

注意該噪聲為matlab中的椒鹽噪聲，可以發現對椒鹽噪聲的處理效果優秀。

最後在此處收錄一份別人寫的matlab程式：

function f = adpmedian(g, Smax)
%ADPMEDIAN Perform adaptive median filtering.
%   F = ADPMEDIAN(G, SMAX) performs adaptive median filtering of
%   image G.  The median filter starts at size 3-by-3 and iterates up
%   to size SMAX-by-SMAX. SMAX must be an odd integer greater than 1.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2003/11/21 14:19:05 $

% SMAX must be an odd, positive integer greater than 1.
if (Smax <= 1) | (Smax/2 == round(Smax/2)) | (Smax ~= round(Smax))
   error('SMAX must be an odd integer > 1.')
end
[M, N] = size(g);

% Initial setup.
f = g;
f(:) = 0;
alreadyProcessed = false(size(g));

% Begin filtering.
for k = 3:2:Smax
   zmin = ordfilt2(g, 1, ones(k, k), 'symmetric');
   zmax = ordfilt2(g, k * k, ones(k, k), 'symmetric');
   zmed = medfilt2(g, [k k], 'symmetric');
   
   processUsingLevelB = (zmed > zmin) & (zmax > zmed) & ...
       ~alreadyProcessed; 
   zB = (g > zmin) & (zmax > g);
   outputZxy  = processUsingLevelB & zB;
   outputZmed = processUsingLevelB & ~zB;
   f(outputZxy) = g(outputZxy);
   f(outputZmed) = zmed(outputZmed);
   
   alreadyProcessed = alreadyProcessed | processUsingLevelB;
   if all(alreadyProcessed(:))
      break;
   end
end

% Output zmed for any remaining unprocessed pixels. Note that this
% zmed was computed using a window of size Smax-by-Smax, which is
% the final value of k in the loop.
f(~alreadyProcessed) = zmed(~alreadyProcessed);
 

參考資源：