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

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

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