影像處理–中值濾波

中值濾波

 中值濾波就是是基於排序統計理論的一種能有效抑制雜訊的非線性訊號處理技術，它的基本原理是把數位影像或數字序列中一點的值用該點的一個鄰域中各點值的中值代替，讓周圍的像素值接近的值，從而消除孤立的雜訊點。中值濾波一般使用模板的方法實現，對模板內的像素按照像素值的大小進行排序，產生單調上升（或下降）的二維資料序列，並使用下面的公式進行輸出：

         g（x,y）=med{f(x-m,y-n),(m,n∈W)}

其中,f(x,y)表示原始的映像，而g(x,y)表示處理後的映像。

        中值濾波一般使用二維模板，濾波視窗通常為 3*3，5*5，7*7 地區，實際使用中，我們常常放大視窗長度，選取最合適的直到濾波效果滿意為止，對於緩變長輪廓物體一般採用方形和圓形，對於尖角形一般採用十字形視窗。後面的程式採用 3*3 矩形地區。實現方法是通過從映像中的某個採樣視窗取出奇數個資料進行排序。用排序後的中值取代要處理的資料即可。

中值濾波適用範圍

中值濾波對椒鹽雜訊的抑制效果好，在抑制隨機雜訊的同時能有效保護邊緣少受模糊。但它對點、線等細節較多的映像卻不太合適。對中值濾波法來說，正確選擇視窗尺寸的大小是很重要的環節。一般很難事先確定最佳的視窗尺寸，需通過從小視窗到大視窗的中值濾波實驗，再從中選取最佳的。

在opencv下使用中值濾波

     OpenCV函數 medianBlur 執行中值濾波操作:

/*************************************************************************   *   *   函數名稱：   *       MedianFilter()   *   *   參數:   *       LPSTR   lpDIBBits -   指向源DIB映像指標   *       LONG     lWidth -   源映像寬度（象素數）   *       LONG     lHeight -   源映像高度（象素數）   *       int       iFilterH -   濾波器的高度   *       int       iFilterW -   濾波器的寬度   *       int       iFilterMX -   濾波器的中心元素X座標   *       int       iFilterMY -   濾波器的中心元素Y座標   *   *   傳回值:   *       BOOL -   成功返回TRUE，否則返回FALSE。   *   *   說明:   *       該函數對DIB映像進行中值濾波。   *   ************************************************************************/   BOOL   WINAPI   MedianFilter(LPSTR   lpDIBBits,   LONG   lWidth,   LONG   lHeight,     int   iFilterH,   int   iFilterW,     int   iFilterMX,   int   iFilterMY)   {   //   指向源映像的指標   unsigned   char* lpSrc;   //   指向要複製地區的指標   unsigned   char* lpDst;   //   指向複製映像的指標   LPSTR lpNewDIBBits;   HLOCAL hNewDIBBits;   //   指向濾波器數組的指標   unsigned   char *   aValue;   HLOCAL hArray;   //   迴圈變數   LONG i;   LONG j;   LONG k;   LONG l;   //   映像每行的位元組數   LONG lLineBytes;   //   計算映像每行的位元組數   lLineBytes   =   WIDTHBYTES(lWidth   *   8);   //   暫時分配記憶體，以儲存新映像   hNewDIBBits   =   LocalAlloc(LHND,   lLineBytes   *   lHeight);   //   判斷是否記憶體配置失敗   if   (hNewDIBBits   ==   NULL)   {   //   分配記憶體失敗   return   FALSE;   }   //   鎖定記憶體   lpNewDIBBits   =   (char   *   )LocalLock(hNewDIBBits);   //   初始化映像為原始映像   memcpy(lpNewDIBBits,   lpDIBBits,   lLineBytes   *   lHeight);   //   暫時分配記憶體，以儲存濾波器數組   hArray   =   LocalAlloc(LHND,   iFilterH   *   iFilterW);   //   判斷是否記憶體配置失敗   if   (hArray   ==   NULL)   {   //   釋放記憶體   LocalUnlock(hNewDIBBits);   LocalFree(hNewDIBBits);   //   分配記憶體失敗   return   FALSE;   }   //   鎖定記憶體   aValue   =   (unsigned   char   *   )LocalLock(hArray);   //   開始中值濾波   //   行(除去邊緣幾行)   for(i   =   iFilterMY;   i   <   lHeight   -   iFilterH   +   iFilterMY   +   1;   i++)   {   //   列(除去邊緣幾列)   for(j   =   iFilterMX;   j   <   lWidth   -   iFilterW   +   iFilterMX   +   1;   j++)   {   //   指向新DIB第i行，第j個象素的指標   lpDst   =   (unsigned   char*)lpNewDIBBits   +   lLineBytes   *   (lHeight   -   1   -   i)   +   j;   //   讀取濾波器數組   for   (k   =   0;   k   <   iFilterH;   k++)   {   for   (l   =   0;   l   <   iFilterW;   l++)   {   //   指向DIB第i   -   iFilterMY   +   k行，第j   -   iFilterMX   +   l個象素的指標   lpSrc   =   (unsigned   char*)lpDIBBits   +   lLineBytes   *   (lHeight   -   1   -   i   +   iFilterMY   -   k)   +   j   -   iFilterMX   +   l;   //   儲存象素值   aValue[k   *   iFilterW   +   l]   =   *lpSrc;   }   }   //   擷取中值   *   lpDst   =   GetMedianNum(aValue,   iFilterH   *   iFilterW);   }   }   //   複製變換後的映像   memcpy(lpDIBBits,   lpNewDIBBits,   lLineBytes   *   lHeight);   //   釋放記憶體   LocalUnlock(hNewDIBBits);   LocalFree(hNewDIBBits);   LocalUnlock(hArray);   LocalFree(hArray);   //   返回   return   TRUE;   }   /*************************************************************************   *   *   函數名稱：   *       GetMedianNum()   *   *   參數:   *       unsigned   char   *   bpArray -   指向要擷取中值的數組指標   *       int       iFilterLen -   數組長度   *   *   傳回值:   *       unsigned   char             -   返回指定數組的中值。   *   *   說明:   *       該函數用冒泡法對一維數組進行排序，並返回數組元素的中值。   *   ************************************************************************/   unsigned   char   WINAPI   GetMedianNum(unsigned   char   *   bArray,   int   iFilterLen)   {   //   迴圈變數   int i;   int j;   //   中間變數   unsigned   char   bTemp;   //   用冒泡法對數組進行排序   for   (j   =   0;   j   <   iFilterLen   -   1;   j   ++)   {   for   (i   =   0;   i   <   iFilterLen   -   j   -   1;   i   ++)   {   if   (bArray[i]   >   bArray[i   +   1])   {   //   互換   bTemp   =   bArray[i];   bArray[i]   =   bArray[i   +   1];   bArray[i   +   1]   =   bTemp;   }   }   }   //   計算中值   if   ((iFilterLen   &   1)   >   0)   {   //   數組有奇數個元素，返回中間一個元素   bTemp   =   bArray[(iFilterLen   +   1)   /   2];   }   else   {   //   數組有偶數個元素，返回中間兩個元素平均值   bTemp   =   (bArray[iFilterLen   /   2]   +   bArray[iFilterLen   /   2   +   1])   /   2;   }   //   返回中值   return   bTemp;   }

CUDA並行條件下的中值濾波

從中值濾波的原理可以看出,平滑後映像的每個像素的值只與原圖該點其鄰域的像素值有關,並且每個像素值的中值濾波處理演算法是完全相同的。並且演算法為只涉及局部運算

的像素級處理，並且處理的資料量巨大，局部資料之間的相關性小，沒有涉及知識推理和人工幹預，演算法本身的並行化程度很高，並且 CUDA 的 SIMT（單指令多線程）並行流式程式處理模式在進行這樣的像素級的影像處理時具有明顯的優勢。設計好核心運行函數，使平滑映像中的每一個像素對應一個由一個線程產生。由此首先設計 kernel 函數MedianFilter 。

_global_ void MedianFilter(int In[N][N],int Out[N][N],int Width,int Height){int window[9];unsigned int x=blockIdx.x*blockDim.x+threadIdx.x;unsigned int y=blockIdx.y*blockDim.y+threadIdx.y;if(x>= Width && y>= Height) return;window[0]=(y==0||x==0)?0:In[(y-1)* Width +x-1];window[1]=(y==0)?0:In[(y-1)* Width +x];window[2]=(y==0||x==DATA_W-1)? 0:In[(y-1)* Width +x+1];window[3]=(x==0)? 0:In[y* Width +x-1];window[4]= In[y* Width +x];window[5]=(x==DATA_W-1)? 0:In[y* Width +x+1];window[6]=(y==DATA_H-1||x==0)? 0:In[(y+1)* Width +x-1];window[7]=(y==DATA_H-1)? 0:In[(y+1)* Width +x];window[8]=(y==DATA_H-1||x==DATA_W-1)? 0:In[(y+1)* Width +x+1];for (unsigned int j=0; j<5; ++j){int min=j;for (unsigned int l=j+1; l<9; ++l)if (window[l] < window[min])min=l;const float temp=window[j];window[j]=window[min];window[min]=temp;}Out[y* Width + x]=window[4];}

說明：

       根據 CUDA 程式定義的文法，_global_ 修飾符聲明函數 MedianFilter 為kernel 函數,表明此類函數在裝置上執行，僅可通過宿主調用。函數 MedianFilter 有四個
形參：int In[N][N]定義了指向輸入的待處理的映像資料指標, int Out[N][N] 定義了指向輸出的處理後的映像資料指標,int Width 定義了待處理的映像的寬度,int Height 定義
了待處理的映像的高度。int window[9]表示為線程定義了一個整形數組用於排序。int x和 y 求出當前線程的座標，確定我們用該線程在產生結果像素的座標。if 的判斷線程超出
邊界說明座標不合法，退出。window[0]到 window[8] 讀取當前點及其 8 鄰域像素值。for迴圈對當前點及其 8 鄰域像素值排序，最後輸出中值濾波結果。

主機端代碼如下

void main ( ){ //讀入待銳利化影像檔並初始化主機;//將映像資料轉送到裝置;dim3 dimBlock (16,16) ;dim3 dimGrid(（Width+dimBlock.x–1）/dimBlock.x,（Height+dimBlock.y–1)/dimBlock.y);clock_t begin = clock ( ) ;, end;MedianFilter<<<dimGrid,dimBlock>>>(inlp,outlp,Width, Height);cudaThreadSynchronize();double cost;end = clock();cost = (double)(end - begin) / CLOCKS_PER_SEC;}