Delphi图像处理 -- 高斯模糊

阅读提示：

    《Delphi图像处理》系列以效率为侧重点，一般代码为PASCAL，核心代码采用BASM。

    《C++图像处理》系列以代码清晰，可读性为主，全部使用C++代码。

    尽可能保持二者内容一致，可相互对照。

    本文代码必须包括文章《Delphi图像处理 -- 数据类型及公用过程》中的ImageData.pas单元。

 

    说明：图像高斯模糊处理代码修改次数最多，此次的修改虽然没有改变算法，但是处理流程做了修改，仅此就可以在原有基础上提高速度40%以上的。同时，这次采用了SSE浮点运算替代了原来一般汇编的定点数运算。为了方便比较，同时也不想毁去原有代码，将原代码略作修改后继续保留，此次修改的代码附在文章后面。

     我在文章《Delphi图像处理 -- 图像卷积》中，曾经介绍过利用通用的图像卷积过程对图像进行高斯模糊处理，其处理效果还不错，处理小型图像时感觉也还行，但是处理较大图像时的速度还是嫌慢，在我的P4 2.8G、1G内存的机器上对千万像素图像进行Q=3，R=5的高斯模糊处理，不包括图像装载和前期数据转换，耗时达8600ms以上，虽经几次修改，其处理速度始终得不到明显提高，主要原因还是采用通用卷积过程处理的问题：用R=5得到的卷积模板为11*11像素，一个像素有4个分量（32位ARGB），对每个象素必须作11*11*4=484个乘法、484个加法及4个除法，最后还得作4个分量是否超界的判断处理，想不慢都难啦！如果不是采用BASM定点数处理代码，其处理速度更是难以想象。

    我在网上多次查找图像高斯模糊的优化算法，不少算法和处理方式，包括代码优化还不如我的那个高斯模糊处理过程，使我很失望。前天查找其它资料时，在国外某个网站上发现介绍图像高斯模糊处理方法时似乎与常规的算法有所不同，但都没有详细的资料（因为不懂外语，很少上国外网站，但看些公式、伪代码还是行的）， 经过反复琢磨，可以将其处理流程归纳如下：

    1、用给定的确定Q和长度size，计算一个size+1长的高斯分布权数数据weights： 

// 计算初始数据for (i = -radius; i <= radius; i ++){    x = i / Q;    weights[i+radius] = exp(-x * x / 2)}// 求和sum = 0for (i = -radius; i <= radius; i ++){    sum += weights[i+radius]}// 数据归一，即归一后的数据之和等于1for (i = -radius; i <= radius; i ++){    weights[i+radius] /= sum}

    2、使用weights对原图像作垂直的模糊运算，即以像素(x, y)为中心，对(x, y - radius)和(x, y + radius)的像素点与weights对应的值相乘后求和得到新的像素，并写入到一个临时的图像上相应的点上（因为数据进行了归一处理，乘积和不必再作除法运算）；

    3、使用weights对临时图像作水平的模糊运算，即以像素(x, y)为中心，对(x - radius, y)和(x + radius, y)的像素点与weights对应的相乘后求和得到新的像素，并写入到目标图像上相应的点上。

    处理过程结束。

    由于上面的处理流程只是对图像每个象素作了一个“十”字型的运算，使得对每个象素点的运算大大减少，模糊长度越大，减少的越多。如前面所说的Q=3、R=5的模糊运算只需要11*2*4=88个乘法、88个加法即可。

    我还是采用BASM按上面的流程作定点数运算，改进后的高斯模糊过程代码如下：

 

procedure CrossBlur(var Dest: TImageData; const Source: TImageData; Weights: Pointer; Size: Integer);var  height, srcStride: Integer;  _weights: Pointer;  dstOffset, srcOffset: Integer;  reds, greens, blues: Integer;asm    push    esi    push    edi    push    ebx    mov     _Weights, ecx    mov     ecx, [edx].TImageData.Stride    mov     srcStride, ecx    call    _SetCopyRegs    mov     height, edx    mov     dstOffset, ebx    push    esi    push    edi    push    edx    push    ecx    push    eax    // blur col    add     ecx, Size           // width = Source.Width    dec     ecx    mov     edi, _weights       // edi = weights@@cyLoop:    push    ecx@@cxLoop:    push    ecx    push    esi    push    edi    xor     ebx, ebx    mov     reds, ebx    mov     greens, ebx    mov     blues, ebx    mov     ecx, Size@@cblurLoop:    movzx   eax, [esi].TARGBQuad.Blue    movzx   edx, [esi].TARGBQuad.Green    imul    eax, [edi]    imul    edx, [edi]    add     blues, eax    add     greens, edx    movzx   eax, [esi].TARGBQuad.Red    movzx   edx, [esi].TARGBQuad.Alpha    imul    eax, [edi]    imul    edx, [edi]    add     reds, eax    add     ebx, edx    add     edi, 4    add     esi, srcStride    loop    @@cblurLoop    pop     edi    pop     esi    mov     eax, blues    mov     edx, greens    mov     ecx, reds    shr     eax, 16    shr     edx, 16    shr     ecx, 16    shr     ebx, 16    mov     [esi].TARGBQuad.Blue, al    mov     [esi].TARGBQuad.Green, dl    mov     [esi].TARGBQuad.Red, cl    mov     [esi].TARGBQuad.Alpha, bl    add     esi, 4    pop     ecx    loop    @@cxLoop    pop     ecx    dec     height    jnz     @@cyLoop    pop     srcOffset    pop     ecx    pop     height    pop     edi    pop     esi    // blur row@@ryLoop:    push    ecx@@rxLoop:    push    ecx    push    esi    push    edi    xor     ebx, ebx    mov     reds, ebx    mov     greens, ebx    mov     blues, ebx    mov     ecx, Size    mov     edi, _weights@@rblurLoop:    movzx   eax, [esi].TARGBQuad.Blue    movzx   edx, [esi].TARGBQuad.Green    imul    eax, [edi]    imul    edx, [edi]    add     blues, eax    add     greens, edx    movzx   eax, [esi].TARGBQuad.Red    movzx   edx, [esi].TARGBQuad.Alpha    imul    eax, [edi]    imul    edx, [edi]    add     reds, eax    add     ebx, edx    add     edi, 4    add     esi, 4    loop    @@rblurLoop    pop     edi    pop     esi    mov     eax, blues    mov     edx, greens    mov     ecx, reds    shr     eax, 16    shr     edx, 16    shr     ecx, 16    shr     ebx, 16    mov     [edi].TARGBQuad.Blue, al    mov     [edi].TARGBQuad.Green, dl    mov     [edi].TARGBQuad.Red, cl    mov     [edi].TARGBQuad.Alpha, bl    add     esi, 4    add     edi, 4    pop     ecx    loop    @@rxLoop    add     esi, srcOffset    add     edi, dstOffset    pop     ecx    dec     height    jnz     @@ryLoop    pop     ebx    pop     edi    pop     esiend;procedure ImageGaussiabBlur(var Data: TImageData; Q: double; Radius: Integer);var  src: TImageData;  fweights: array of Single;  weights: array of Integer;  i, size: Integer;  fx: Double;begin  if Radius <= 0 then  begin    if Abs(Q) < 1.0 then Radius := 1    else Radius := Round(Abs(Q)) + 2;  end;  size := Radius shl 1 + 1;  SetLength(fweights, size);  for i := 1 to Radius do  begin    fx := i / Q;    fweights[Radius + i] := exp(-fx * fx / 2);    fweights[Radius - i] := fweights[Radius + i];  end;  fweights[Radius] := 1.0;  fx := 0.0;  for i := 0 to size - 1 do    fx := fx + fweights[i];  SetLength(weights, size);  for i := 0 to size - 1 do    weights[i] := Round(fweights[i] / fx * 65536.0);  SetLength(fweights, 0);  src := _GetExpandData(Data, Radius);  CrossBlur(Data, src, weights, size);  FreeImageData(src);end;

    用改进后的高斯模糊处理过程在我的机器上对千万像素图像进行Q=3，R=5的高斯模糊处理，不包括图像装载和前期数据转换，耗时为1390ms，处理速度确实得到了大幅度的提高。我是按32位ARGB颜色处理图像像素的，如果改为24位RGB颜色处理图像像素，耗时还可以减少，不过，RGB颜色没法处理PNG等32位像素格式的图像。

    不用模板卷积方式，而采用“十”字运算进行高斯模糊处理，效果如何呢？请看下面的简单例子代码及处理效果图：

    例子代码： 

procedure TForm1.Button3Click(Sender: TObject);var  bmp: TGpBitmap;  g: TGpGraphics;  data: TImageData;begin  bmp := TGpBitmap.Create('..\media\56-3.jpg');  g := TGpGraphics.Create(Canvas.Handle);  g.DrawImage(bmp, 0, 0);  data := LockGpBitmap(bmp);  ImageGaussiabBlur(Data, 3, 6);  UnlockGpBitmap(bmp, data);  g.DrawImage(bmp, data.Width, 0);  g.Free;  bmp.Free;end;

    处理原图：

    处理效果与Photoshop高斯模糊处理对比图：

    左上是Photoshop半径3.0高斯模糊效果图，右上是本文过程Q=3.0,R=6高斯模糊效果图。

    左下是Photoshop半径5.0高斯模糊效果图，右下是本文过程Q=5.0,R=9高斯模糊效果图。

    怎么样，效果还不错吧！

    遗憾的是我没能找到按照Q自动计算模糊半径的方法，所以处理过程给出了2个参数Q和Radius。

    下面是本次修改后的SSE代码，因原理和算法同上，只是在处理手法上有些不同：因为高斯模糊矩阵上下、左右都是对称的，因此以半径点位中心，将上下对称行（列处理时）或者左右对称列（行处理时）相加后再与高斯分布权数数据相乘，如此，除中心行（列）外，只须作以前的50%处理。

procedure CrossBlur(var Dest: TImageData; const Source: TImageData; Weights: Pointer; Radius: Integer);var  height, srcStride: Integer;  dstOffset, srcOffset: Integer;asm    push      esi    push      edi    push      ebx    push      ecx    mov       ecx, [edx].TImageData.Stride    mov       srcStride, ecx    call      _SetCopyRegs    mov       height, edx    mov       srcOffset, eax    mov       dstOffset, ebx    pop       ebx    pxor      xmm7, xmm7    push      esi           // pst = Source.Scan0    push      edi    push      edx    push      ecx    // blur col    mov       eax, srcStride    mov       edx, eax    shr       edx, 2        // width = Source.Width    mov       edi, Radius    shl       edi, 1    imul      edi, eax    add       edi, esi      // psb = pst + Radius * 2 * Source.Stride@@cyLoop:    push      edx@@cxLoop:    push      esi    push      edi    push      ebx    mov       ecx, Radius    pxor      xmm0, xmm0    // sum = 0@@cblurLoop:    movd      xmm1, [esi]   // for (i = 0; i < Radius; i ++)    movd      xmm2, [edi]   // {    punpcklbw xmm1, xmm7    punpcklbw xmm2, xmm7    paddw     xmm1, xmm2    //   ps = pst + psb    punpcklwd xmm1, xmm7    cvtdq2ps  xmm1, xmm1    //   pfs (flaot * 4) = ps (int * 4)    mulps     xmm1, [ebx]   //   pfs *= Weights[i]    addps     xmm0, xmm1    //   sum += pfs    add       ebx, 16    add       esi, eax      //   pst += Source.Stride    sub       edi, eax      //   psb -= Source.Stride    loop      @@cblurLoop   // }    movd      xmm1, [esi]    punpcklbw xmm1, xmm7    punpcklwd xmm1, xmm7    cvtdq2ps  xmm1, xmm1    // pfs (flaot * 4) = pst (int * 4)    mulps     xmm1, [ebx]   // pfs *= Weights[Radius]    addps     xmm0, xmm1    // sum += pfs    pop       ebx    pop       edi    pop       esi    cvtps2dq  xmm0, xmm0    // ps (int * 4) = sum (flaot * 4)    packssdw  xmm0, xmm7    packuswb  xmm0, xmm7    movd      [esi], xmm0   // pst (byte * 4) = ps (int * 4) pask    add       esi, 4    add       edi, 4    dec       edx    jnz       @@cxLoop    pop       edx    dec       height    jnz       @@cyLoop    pop       edx    pop       height    pop       edi           // pd = Dest.Scan0    pop       esi           // psl = pst    mov       eax, Radius    shl       eax, 1+2    add       eax, esi      // psr = psl + Radius * 2    // blur row@@ryLoop:    push      edx           // width = Dest.Width@@rxLoop:    push      esi    push      ebx    push      eax    mov       ecx, Radius    pxor      xmm0, xmm0    // sum = 0@@rblurLoop:    movd      xmm1, [esi]   // for (i = 0; i < Radius; i ++)    movd      xmm2, [eax]   // {    punpcklbw xmm1, xmm7    punpcklbw xmm2, xmm7    paddw     xmm1, xmm2    //   ps = psl + psr    punpcklwd xmm1, xmm7    cvtdq2ps  xmm1, xmm1    //   pfs (flaot * 4) = ps (int * 4)    mulps     xmm1, [ebx]   //   pfs *= Weights[i]    addps     xmm0, xmm1    //   sum += pfs    add       ebx, 16    add       esi, 4        //   psl ++    sub       eax, 4        //   psr --    loop      @@rblurLoop   // }    movd      xmm1, [esi]    punpcklbw xmm1, xmm7    punpcklwd xmm1, xmm7    cvtdq2ps  xmm1, xmm1    // pfs (flaot * 4) = psl (int * 4)    mulps     xmm1, [ebx]   // pfs *= Weights[Radius]    addps     xmm0, xmm1    // sum += pfs    cvtps2dq  xmm0, xmm0    // ps (int * 4) = sum (flaot * 4)    packssdw  xmm0, xmm7    packuswb  xmm0, xmm7    movd      [edi], xmm0   // pd (byte * 4) = ps (int * 4) pask    pop       eax    pop       ebx    pop       esi    add       eax, 4    add       esi, 4    add       edi, 4    dec       edx    jnz       @@rxLoop    add       eax, srcOffset    add       esi, srcOffset    add       edi, dstOffset    pop       edx    dec       height    jnz       @@ryLoop    pop       ebx    pop       edi    pop       esiend;// --> st x// <-- st e**x = 2**(x*log2(e))function _Expon: Extended;asm    fldl2e              // y = x*log2e    fmul    fld     st(0)       // i = round(y)    frndint    fsub    st(1), st   // f = y - i    fxch    st(1)       // z = 2**f    f2xm1    fld1    fadd    fscale              // result = z * 2**i    fstp    st(1)end;function GetWeights(var Buffer, Weights: Pointer; Q: Single; Radius: Integer): Integer;const  _fcd1: Single = 0.1;  _fc1: Single = 1.0;  _fc2: Single = 2.0;  _fc250: Single = 250.0;  _fc255: Single = 255.0;var  R: Integer;  v, QQ2: double;asm    mov     R, ecx    mov     ecx, eax    fld     Q    fabs    fcom    _fcd1    fstsw   ax    sahf    jae     @@1    fld     _fcd1    fstp    st(1)               // if (Q < 0.1) Q = 0.1    jmp     @@2@@1:    fcom    _fc250    fstsw   ax    sahf    jbe     @@2    fld     _fc250    fstp    st(1)               // if (Q > 250) Q = 250@@2:    fst     Q    fmul    Q    fmul    _fc2    fstp    QQ2                 // QQ2 = 2 * Q * Q    fwait    mov     eax, R    test    eax, eax    jg      @@10    push    eax                 // if (radius <= 0)    fld1                        // {    fadd    Q                   //   radius = Abs(Q) + 1    fistp   [esp].Integer    fwait    pop     eax@@testRadius:                   //   while (TRUE)    mov     R, eax              //   {    fldz                        //     sum = 0@@testLoop:                     //     for (R = radius; R > 0; R ++)    fild    R                   //     {    fld     st(0)    fmulp   st(1), st    fdiv    QQ2    fchs    call    _Expon              //       tmp = Exp(-(R * R) / (2.0 * Q * Q));    cmp     R, eax    jne     @@3    fst     v                   //       if (R == radius) v = tmp@@3:    faddp   st(1), st(0)        //       sum += tmp    dec     R    jnz     @@testLoop          //     }    fmul    _fc2                //     sum *= 2    fadd    _fc1                //     sum += 1    fdivr   v    fmul    _fc255    fistp   R    cmp     R, 0    je      @@4                 //     if ((INT)(v / sum * 255 + 0.5) = 0) break    inc     eax                 //     radius ++    jmp     @@testRadius        //   }@@4:    dec     eax    jnz     @@5    inc     eax@@5:    mov     R, eax              // }@@10:    inc     eax    shl     eax, 4    add     eax, 12    push    edx    push    ecx    mov     edx, eax    mov     eax, GHND    call    GlobalAllocPtr    pop     ecx    pop     edx    test    eax, eax    jz      @@Exit    mov     [ecx], eax          // buffer = GlobalAllocPtr(GHND, (Radius + 1) * 16 + 12)    add     eax, 12    and     eax, -16    mov     [edx], eax          // weights = ((char* )buffer + 12) & 0xfffffff0    mov     ecx, R              // ecx = radius    mov     edx, eax            // edx = weights    fldz                        // for (i = radius, sum = 0; i > 0; i --)@@clacLoop:                     // {    fild    R    fld     st(0)    fmulp   st(1), st    fdiv    QQ2    fchs    call    _Expon    fstp    [edx].Double        //   weights[i] = Expon(-(i * i) / (2 * Q * Q))    fadd    [edx].Double        //   sum += weights[i]    add     edx, 16    dec     R    jnz     @@clacLoop          // }    fmul    _fc2                // sum *= 2    fld1    fstp    [edx].Double        // weights[radius] = 1    fadd    [edx].Double        // sum += weights[radius]    push    ecx    inc     ecx@@divLoop:                      // for (i = 0; i <= Radius; i ++)    fld     st(0)               //   weights[i] = Round(weights[i] / sum)    fdivr   [eax].Double    fst     [eax].Single    fst     [eax+4].Single    fst     [eax+8].Single    fstp    [eax+12].Single    add     eax, 16    loop    @@divLoop    ffree   st(0)    fwait    pop     eax                 // return Radius@@Exit:end;procedure ImageGaussiabBlur(var Data: TImageData; Q: Single; Radius: Integer);var  Buffer, Weights: Pointer;  src: TImageData;begin  Radius := GetWeights(Buffer, Weights, Q, Radius);  if Radius = 0 then Exit;  if Data.AlphaFlag then    ArgbConvertPArgb(Data);  src := _GetExpandData(Data, Radius);  CrossBlur(Data, src, Weights, Radius);  FreeImageData(src);  GlobalFreePtr(Buffer);  if Data.AlphaFlag then    PArgbConvertArgb(Data);end;

    《Delphi图像处理》系列使用GDI+单元下载地址和说明见文章《GDI+ for VCL基础 -- GDI+ 与 VCL》。

    因水平有限，错误在所难免，欢迎指正和指导。邮箱地址：[e-mail]

    这里可访问《Delphi图像处理 -- 文章索引》。