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 = 0
for (i = -radius; i <= radius; i ++)
{
    sum += weights[i+radius]
}

// 数据归一，即归一后的数据之和等于1
for (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     esi
end;

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

// --> 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》。

    因水平有限，错误在所难免，欢迎指正和指导。邮箱地址：maozefa@hotmail.com

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