Delphi影像處理 — 中值濾波（灰階分組統計法）

閱讀提示：

    《Delphi影像處理》系列以效率為側重點，一般代碼為PASCAL，核心代碼採用BASM。

    《C++影像處理》系列以代碼清晰，可讀性為主，全部使用C++代碼。

    儘可能保持二者內容一致，可相互對照。

    本文代碼必須包括文章《Delphi影像處理 -- 資料類型及公用過程》中的ImageData.pas單元。

 

    這是《Delphi影像處理 -- 中值濾波》一文的改進版。亦可參見《C++影像處理 -- 中值濾波》。

    中值濾波是影像處理中常用的一種雜訊濾波方法。傳統的映像中值濾波代碼採用排序方法實現，處理速度主要取決於排序演算法，但無論什麼排序演算法，總離不開大量的元素比較、交換或移動，而這些恰好是當前電腦處理的“弱項”（有經驗的程式員都知道，電腦資料處理中，比較、轉移、交換和頻繁的資料移動比直接的算術運算和邏輯運算耗時多了），再加上沒有一種好的排序演算法能同時適應不同濾波半徑的資料排序速度，所以在傳統中值濾波實現代碼中多使用選擇排序、冒泡排序或者直接排序等簡單排序演算法，進階點的如快速排序法用在中值濾波代碼中往往會使處理速度更慢。對於半徑為1的中值濾波倒是有一種較好的排序演算法，我在《Delphi影像處理 -- 中值濾波》一文中實現過，處理速度還是較快的。

    既然排序過程是映像中值濾波處理的瓶頸，能不能拋開它，用其它手段實現呢？這就是本文要探討的問題。有朋友可能會有疑問，不排序怎麼擷取中間值呢，是否採用網上有些文章介紹的近似值來代替？不，本文介紹的方法決不是近似中間值，而是的的確確的“精確”中間值。

    映像中值濾波中的中間值。在統計學中叫做中位元，是平均數指標的一種。平均數指標按統計的複雜程度可分為簡單平均數和加權平均數，所謂簡單平均數就是對統計總體的每個個體進行累計計算後求得；而加權平均數則是先對統計總體的所有個體進行分組，然後以各個組的大小作為權數後進行累計計算所得。中位元既然是平均數指標的一種，當然也存在兩種計算方法。加權中位元和加權算術平均數的區別在於後者是對各個分組用權數相乘進行累積後除以總體個數，而前者只需要對各組權數進行累積後除以2，累積權數大於或等於這個數的組的值就是中位元。傳統中值濾波實現代碼中採用的排序手段實質就是簡單中位平均數統計方法，而本文要介紹的是採用分組加權的方法來擷取中位平均數，即中值濾波的中間值。

    採用分組加權統計方法的前提條件是對無限統計總體進行有限的分組，例如，人口年齡統計，就是首先確定各個年齡段，這個年齡段可以是一歲、五歲、十歲等。而映像像素的R、G、B值在0 -- 255之間，正好符合有限分組的前提條件。說到這裡，很多人可能已經明白我要表達的意思了：

    1、按R、G、B分別定義一個256大小的灰階統計資料；

    2、將映像像素總體的每個個體像素的R、G、B值進行歸類；

    3、確定累積中間值權數的大小；

    4、從數組左端或者右端開始權數累計，累積到大於或等於中間值權數的那個灰階組就是中間值！

    從前面幾個步驟不難看出，前3個步驟就是影像處理中灰階統計的分類方法，第4個累積步驟與灰階統計累積計算有2個不同點，一是灰階統計累積的是權數*灰階，而這裡只累積灰階；二是灰階統計累積需要全部完成，而中值累積只要達到中間值權數那個組就可以終止。在這4個步驟中，前3個步驟是相當快的（實際上只是第二個步驟），因為其中既無乘除運算，也沒有比較轉移，更沒有元素交換或移動等動作。而制約資料處理速度的瓶頸就在第4步，因為對每個像素都必須分別按R、G、B通道對一共768個元素大小的資料進行累積，哪怕是每個通道除了2個加法運算和唯一的一次比較判斷外，沒有其它運算，也不需要每次都必須累積到位，但仍然是比較耗時的，不過本文在代碼實現過程中，儘可能地作了一些彌補。最後結果同排序方法比起來，可算是相當快捷了，而且濾波半徑越大，差距越明顯，幾十倍的差距絕不是天方夜譚！就算是同我前面所說的半徑為1的改進排序演算法比較來，大多數情況下，也略有勝出。

procedure MedianValue(var Dest: TImageData; const Source: TImageData;  MedianGray, Size, Stride: Integer);var  buffer: array[0..767] of Integer;  redAddr, greenAddr, blueAddr, delta: Integer;  redOff, greenOff, blueOff: Integer;  width, height, dstOffset, srcOffset: Integer;  median, rowOffset, sizeOffset, mOffset: Integer;asm    push    esi    push    edi    push    ebx    push    ecx    call    _SetCopyRegs    mov     width, ecx    mov     height, edx    mov     dstOffset, ebx    add     eax, 4    mov     srcOffset, eax    pop     ecx    mov     eax, Size    mov     edx, Stride    mov     ebx, edx    imul    edx, eax    shl     eax, 2    mov     sizeOffset, eax    // sizeOffset = Size * 4    sub     ebx, eax    mov     rowOffset, ebx     // rowOffset = Stride - Size * 4    sub     edx, 4    mov     mOffset, edx       // mOffset = Stride * Size - 4    mov     eax, Size    mul     Size    inc     eax    shr     eax, 1    mov     median, eax         // median = (size * size + 1) / 2    mov     blueAddr, 0    mov     greenAddr, 256*4    mov     redAddr, 512*4    cmp     ecx, 128    jb      @@1    mov     blueOff, 256*4    mov     greenOff, 512*4    mov     redOff, 768*4    mov     delta, -4    jmp     @@2@@1:    mov     blueOff, -4    mov     greenOff, 256*4-4    mov     redOff, 512*4-4    mov     delta, 4@@2:    lea     eax, buffer    add     blueAddr, eax    add     greenAddr, eax    add     redAddr, eax    add     blueOff, eax    add     greenOff, eax    add     redOff, eax@@yLoop:    push    width    // RGB灰階統計緩衝區清零    push    edi    mov     edi, blueAddr    mov     ecx, 768    xor     eax, eax    rep     stosd    pop     edi    // 對每行第一個像素臨近地區(Size*Size)的RGB進行灰階統計    push    esi    mov     ecx, Size@@statY:    push    ecx    mov     ecx, Size@@statX:    movzx   eax, [esi].TARGBQuad.Blue    movzx   edx, [esi].TARGBQuad.Green    movzx   ebx, [esi].TARGBQuad.Red    inc     dword ptr buffer[eax*4]    inc     dword ptr buffer[edx*4+256*4]    inc     dword ptr buffer[ebx*4+512*4]    add     esi, 4    loop    @@statX    pop     ecx    add     esi, rowOffset    loop    @@statY    pop     esi    jmp     @@setValue@@xLoop:    // 剔除上個座標點最左邊一列像素的RGB的統計值，    // 追加當前像素最右邊一列像素RGB的統計值    mov     ecx, Size    mov     ebx, sizeOffset    add     ebx, esi@@subLoop:    movzx   eax, [esi].TARGBQuad.Blue    movzx   edx, [esi].TARGBQuad.Green    dec     dword ptr buffer[eax*4]    dec     dword ptr buffer[edx*4+256*4]    movzx   eax, [esi].TARGBQuad.Red    movzx   edx, [ebx].TARGBQuad.Blue    dec     dword ptr buffer[eax*4+512*4]    inc     dword ptr buffer[edx*4]    movzx   eax, [ebx].TARGBQuad.Green    movzx   edx, [ebx].TARGBQuad.Red    inc     dword ptr buffer[eax*4+256*4]    inc     dword ptr buffer[edx*4+512*4]    add     esi, Stride    add     ebx, Stride    loop    @@subLoop    sub     esi, mOffset@@setValue:    mov     edx, median    mov     ecx, delta    mov     ebx, blueOff    xor     eax, eax@@blueLoop:    add     ebx, ecx    add     eax, [ebx]    cmp     eax, edx    jb      @@blueLoop    sub     ebx, blueAddr    shr     ebx, 2    mov     [edi].TARGBQuad.Blue, bl    mov     ebx, greenOff    xor     eax, eax@@greenLoop:    add     ebx, ecx    add     eax, [ebx]    cmp     eax, edx    jb      @@greenLoop    sub     ebx, greenAddr    shr     ebx, 2    mov     [edi].TARGBQuad.Green, bl    mov     ebx, redOff    xor     eax, eax@@redLoop:    add     ebx, ecx    add     eax, [ebx]    cmp     eax, edx    jb      @@redLoop    sub     ebx, redAddr    shr     ebx, 2    mov     [edi].TARGBQuad.Red, bl    add     edi, 4    dec     width    jnz     @@xLoop@@xEnd:    add     esi, srcOffset    add     edi, dstOffset    pop     width    dec     height    jnz     @@yLoop    pop     ebx    pop     edi    pop     esiend;function GetMedianGray(const Data: TImageData): Integer;var  buffer: array[0..255] of Integer;asm    push    esi    push    edi    push    ebx    mov     esi, eax    lea     edi, buffer    mov     ecx, 256    xor     eax, eax    rep     stosd    mov     ecx, [esi].TImageData.Width    mov     edx, [esi].TImageData.Height    mov     ebx, [esi].TImageData.Stride    mov     edi, [esi].TImageData.Scan0    shr     ecx, 3    shr     edx, 3    sal     ebx, 3    mov     eax, ecx    shl     eax, 3+2    sub     ebx, eax    lea     esi, buffer    push    ebp    mov     ebp, edx    mov     eax, ecx    mul     ebp    lea     eax, [eax+eax*2]    shr     eax, 1            // MedianValue = data.width * data.height * 3 / 2    push    eax@@yLoop:    push    ecx@@xLoop:    movzx   eax, [edi].TARGBQuad.Blue    movzx   edx, [edi].TARGBQuad.Green    inc     dword ptr [esi+eax*4]    inc     dword ptr [esi+edx*4]    movzx   eax, [edi].TARGBQuad.Red    inc     dword ptr [esi+eax*4]    add     edi, 32    dec     ecx    jg      @@xLoop    pop     ecx    add     edi, ebx    dec     ebp    jg      @@yLoop    pop     edx    pop     ebp    mov     edi, esi    add     edi, 256*4    mov     ecx, -4    xor     eax, eax@@stat:    add     edi, ecx    add     eax, [edi]    cmp     eax, edx    jb      @@stat    mov     eax, edi    sub     eax, esi    shr     eax, 2    pop     ebx    pop     edi    pop     esiend;procedure ImageMedianValue(var Data: TImageData; Radius: Integer);var  exp: TImageData;begin  if Radius <= 0 then Radius := 1;  exp := _GetExpandData(Data, Radius);  MedianValue(Data, exp, GetMedianGray(Data), (Radius shl 1) + 1, exp.Stride);  FreeImageData(exp);end;

    實現代碼中，在前面所說的4個步驟基礎上作了3點完善：

    1、並非對映像的所有像素都按4個步驟進行。除了每行行首像素作了完整的分類統計外，其它像素只是從統計數組中剔除前一像素臨近值的最左邊一列和追加當前像素臨近值的最右邊一列，這無疑加快了分類統計速度，而且濾波半徑越大，效果越明顯。

    2、在對映像進行濾波處理前，對整個映像像素總體的1/8樣本求了一次中間值權數，其作用是確定像素分類後的累積方向，如果這個總體中間值權數小於128，從灰階統計資料左邊（低端）開始累積，否則則從灰階統計資料右邊（高端）開始累積。通過這個總體中間值權數就可以大致確定該影像處理速度的快慢，如果這個值接近兩端，處理速度相對就快些，反之如靠近128附近，處理速度則相應會慢些，同樣大小的圖片，因這個原因，可能造成成倍的差距，但是最壞的情況下，對灰階統計數組累積時的平均值也不會超過一半。

    3、製作濾波備份源圖時按濾波半徑對映像邊緣作了擴充，這樣有利於對映像邊緣進行濾波處理。有些人在寫中值濾波代碼時往往對映像邊緣略過不作處理，這在小半徑濾波範圍內還無所謂，但是對於較大半徑的濾波處理後，由於未進行濾波處理的邊緣較寬。會使得映像很難看。

    正式由於第2點中所說的差距原因，本文也不好給出具體的測試資料，但有一點是肯定的：同等條件下（語言、程式員水平及編譯器最佳化程度等），比傳統排序法實現的中值濾波要快不少，濾波半徑越大，差距越大！

    最後還是貼個中值濾波測試代碼和：前面是源圖，後面是50半徑中值濾波圖：

procedure TForm1.Button3Click(Sender: TObject);var  bmp: TGpBitmap;  g: TGpGraphics;  data: TImageData;begin  bmp := TGpBitmap.Create('..\media\source.jpg');  g := TGpGraphics.Create(Canvas.Handle);  g.DrawImage(bmp, 0, 0);  data := LockGpBitmap(bmp);  ImageMedianValue(data, 50);  UnlockGpBitmap(bmp, data);  g.DrawImage(bmp, data.Width, 0);  g.Free;  bmp.Free;end;

 

    《Delphi影像處理》系列使用GDI+單元和說明見文章《GDI+ for VCL基礎 -- GDI+ 與 VCL》。

    因水平有限，錯誤在所難免，歡迎指正和指導。郵箱地址：maozefa@hotmail.com

    這裡可訪問《Delphi影像處理 -- 文章索引》。