Delphi影像處理 — USM銳利化

閱讀提示：

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

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

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

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

 

    USM銳利化是用來銳利化映像邊緣的，它通過調整映像邊緣細節的對比，並在邊緣的兩側產生一條亮線一條暗線，使畫面整體更加清晰。

    USM銳利化用公式描述很麻煩，這裡乾脆實現步驟列於下面：

    1、備份映像原資料；

    2、按給定半徑對映像進行高斯模糊；

    3、用原像素與高斯模糊後的像素相減，形成一個差值；

    4、將差值乘以數量；

    5、將差值分為正負2部分，取絕對值；

    6、正負差值分別減去給定的閾值；

    7、原像素加上正差值減去負差值，銳利化完畢。

    下面是USM銳利化代碼，包括高斯模糊代碼（關於高斯模糊見文章《Delphi影像處理 -- 高斯模糊》，下面的高斯模糊代碼也是從該文拷貝來的）：

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;procedure DoUSMSharpen(var Dest: TImageData; const Source: TImageData;  Amount512, Threshold: Integer);asm    push      esi    push      edi    push      ebx    pxor      mm7, mm7    movd      mm6, Threshold    movd      mm5, ecx    pshufw    mm6, mm6, 11000000b    pshufw    mm5, mm5, 11000000b    call      _SetCopyRegs@@yLoop:    push      ecx@@xLoop:    movd      mm1, [esi]  // mm0 = mm1 = 原像素    movd      mm2, [edi]  // mm2 = 高斯模糊像素    movq      mm0, mm1    punpcklbw mm1, mm7    punpcklbw mm2, mm7    psubw     mm1, mm2    psllw     mm1, 7    pmulhw    mm1, mm5    // mm1 = (原像素 - 高斯模糊像素) * 數量    pxor      mm2, mm2    psubw     mm2, mm1    // mm2 = -mm1    psubw     mm1, mm6    // mm1 -= 閾值    psubw     mm2, mm6    // mm2 -= 閾值    packuswb  mm1, mm7    // mm1 = 正離差(小於0的值飽和為0)    packuswb  mm2, mm7    // mm2 = 負離差(小於0的值飽和為0)    paddusb   mm0, mm1    psubusb   mm0, mm2    // mm0 = 原像素 + 正離差 - 負離差    movd      [edi], mm0    add       esi, 4    add       edi, 4    loop      @@xLoop    pop       ecx    add       esi, eax    add       edi, ebx    dec       edx    jnz       @@yLoop    pop       ebx    pop       edi    pop       esi    emmsend;// USM銳利化。參數：映像資料，半徑，數量，閾值procedure ImageUSMSharpen(var Data: TImageData;  Radius: Single; Amount: Integer = 50; Threshold: Integer = 0);  procedure CopyData(var Dest: TImageData; const Source: TImageData);  asm    push    esi    push    edi    mov     ecx, [eax].TImageData.Width    imul    ecx, [eax].TImageData.Height    mov     edi, [eax].TImageData.Scan0    mov     esi, [edx].TImageData.Scan0    rep     movsd    pop     edi    pop     esi  end;var  IsInvertScan0: Boolean;  Src: TImageData;begin  if (Amount < 1) or (Amount > 500) or (not (Threshold in [0..255])) then    Exit;  // 備份映像資料  Src := NewImageData(Data.Width, Data.Height);  IsInvertScan0 := Data.Stride < 0;  if IsInvertScan0 then    _InvertScan0(Data);  CopyData(Src, Data);  // 映像資料作高斯模糊  ImageGaussiabBlur(Data, Radius, 0);  // 用備份映像資料和高斯模糊映像資料作USM銳利化  DoUSMSharpen(Data, Src, (Amount shl 9) div 100, Threshold);  FreeImageData(Src);  if IsInvertScan0 then    _InvertScan0(Data);end;

    下面是個簡單的調用例子：

procedure TForm1.Button2Click(Sender: TObject);var  bmp, bmp2: TBitmap;  data, data2: TImageData;begin  // 從檔案裝入映像資料到bmp，並綁定到data  bmp := TBitmap.Create;  bmp.LoadFromFile('d:\source.bmp');  data := GetBitmapData(bmp);  // USM銳利化：半徑3.0，數量150，閾值10  ImageUSMSharpen(data, 3, 150, 10);  // 畫USM銳利化後的映像  Canvas.Draw(0, 0, bmp);  // 從檔案裝入Photoshop USM銳利化映像資料到bmp2，並綁定到data2  bmp2 := TBitmap.Create;  bmp2.LoadFromFile('d:\sourceUSM.bmp');  data2 := GetBitmapData(bmp2);  // 將ImageUSMSharpen過程銳利化與Photoshop USM銳利化映像資料作比較  ImageCompare(data, data2);  bmp2.Free;  bmp.Free;end;

    下面是原始映像與本文USM銳利化例子運行介面：

    從運行介面上看，本文USM銳利化與Photoshop的USM銳利化有一定的誤差，這個誤差並不大，完全可以接受。但需要說明的是，這裡的誤差是高斯模糊和USM銳利化2個處理結果共同產生的誤差，而《Delphi影像處理 -- 高斯模糊》中高斯模糊實現與Photoshop高斯模糊只是近似，因此，有必要排除高斯模糊影響，來具體確定本文USM銳利化演算法與Photoshop的真實誤差，下面的例子是用Photoshop處理的高斯模糊圖片與原圖作USM銳利化，這時的誤差將完全是本文USM銳利化代碼產生的：

procedure TForm1.Button1Click(Sender: TObject);var  bmp, bmp2: TBitmap;  data, data2: TImageData;begin  // 從檔案裝入Photoshop高斯模糊（半徑3.0）映像資料到bmp，並綁定到data  bmp := TBitmap.Create;  bmp.LoadFromFile('d:\sourceBlur.bmp');  data := GetBitmapData(bmp);  // 從檔案裝入映像資料到bmp2，並綁定到data2  bmp2 := TBitmap.Create;  bmp2.LoadFromFile('d:\source.bmp');  data2 := GetBitmapData(bmp2);  // 用Photoshop高斯模糊映像資料與原資料作USM銳利化：數量150，閾值10  DoUSMSharpen(data, data2, (150 shl 9) div 100, 10);  // 畫USM銳利化後的映像  Canvas.Draw(0, 0, bmp);  // 從檔案裝入Photoshop USM銳利化映像資料到bmp2，並綁定到data2  bmp2.LoadFromFile('d:\sourceUSM.bmp');  data2 := GetBitmapData(bmp2);  // 將DoUSMSharpen過程銳利化與Photoshop USM銳利化映像資料作比較  ImageCompare(data, data2);  bmp2.Free;  bmp.Free;end;

    運行介面：

    運行結果，完全沒有誤差！證明本文USM銳利化演算法和步驟是正確的。

    前面的2個例子中，都調用了映像資料比較過程ImageCompare，也將它列在下面：

procedure ImageCompare(data1, data2: TImageData);var  count, r_count, g_count, b_count: Integer;  diff, r_diff, g_diff, b_diff: Integer;  p1, p2: PARGBQuad;  x, y, offset: Integer;  s: string;begin  r_count := 0;  g_count := 0;  b_count := 0;  r_diff := 0;  g_diff := 0;  b_diff := 0;  p1 := data1.Scan0;  p2 := data2.Scan0;  offset := data1.Stride - data1.Width * sizeof(TARGBQuad);  for y := 1 to data1.Height do  begin    for x := 1 to data1.Width do    begin      diff := p1^.Red - p2^.Red;      if diff <> 0 then      beginInc(r_count);        if diff < 0 then diff := -diff;if r_diff < diff then r_diff := diff;      end;      diff := p1^.Green - p2^.Green;      if diff <> 0 then      beginInc(g_count);if diff < 0 then diff := -diff;if g_diff < diff then g_diff := diff;      end;      diff := p1^.Blue - p2^.Blue;      if diff <> 0 then      beginInc(b_count);if diff < 0 then diff := -diff;if b_diff < diff then b_diff := diff;      end;      Inc(p1);      Inc(P2);    end;    Inc(Integer(p1), offset);    Inc(Integer(p2), offset);  end;  count := data1.Width * data1.Height;  s := Format('像素總數：%d'+ #13 + #10 +              '紅誤差數：%d，誤差率：%d%%，最大誤差：%d' + #13 + #10 +              '綠誤差數：%d，誤差率：%d%%，最大誤差：%d' + #13 + #10 +              '藍誤差數：%d，誤差率：%d%%，最大誤差：%d',              [count, r_count, (r_count * 100) div count, r_diff,               g_count, (g_count * 100) div count, g_diff,               b_count, (b_count * 100) div count, b_diff]);  ShowMessage(s);end;

 

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

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

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