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     esi
end;

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     esi
end;

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图像处理 -- 文章索引》。