Delphi图像处理 -- 平面几何变换（上）

阅读提示：

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

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

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

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

 

    在《Delphi图像处理 -- 平面几何变换类》一文中，介绍了图像平面几何变换类TTransformMatrix，并写了一个简单的临近插值法图像几何变换函数Transform，用于测试。很显然，Transform函数产生的变换图像不仅质量较差，而且也不具备通用性，只能作为一个实现图像几何变换的框架。

    本文拟采用临近插值法、双线性插值法和双立方插值法等三种插值方式，来实现较完整、通用的图形图像平面几何变换。三种插值过程代码包含在《Delphi图像处理 -- 平面几何变换类》一文中，本文代码不包括TTransformMatrix类和三种插值过程。

procedure SetMixerMM;
asm
    pxor      mm7, mm7
    mov       eax, 1011h
    movd      mm6, eax
    pshufw    mm6, mm6, 0
    mov       eax, 8
    movd      mm5, eax
    pshufw    mm5, mm5, 0
end;

procedure MixerColor;
asm
    cmp       eax, 255
    jb        @@ArgbMix
    movd      [edi], xmm0
    ret
@@ArgbMix:
    movdq2q   mm0, xmm0
    movd      mm1, [edi]
    punpcklbw mm0, mm7
    punpcklbw mm1, mm7
    pshufw    mm2, mm0, 255
    movzx     eax, [edi].TARGBQuad.Alpha
    cmp       eax, 255
    jne       @@PArgbMix

    // dest.argb = dest.argb + source.argb - dest.argb * source.alpha / 255
    paddw     mm0, mm1
    pmullw    mm1, mm2
    pmulhuw   mm1, mm6
    paddusw   mm1, mm5
    psrlw     mm1, 4
    psubsw    mm0, mm1
    packuswb  mm0, mm0
    movd      [edi], mm0
    ret
@@PArgbMix:
    // dest.rgb = dest.rgb * dest.alpha / 255
    pshufw    mm3, mm1, 255
    pmullw    mm1, mm3
    pmulhuw   mm1, mm6
    paddusw   mm1, mm5
    psrlw     mm1, 4
    pinsrw    mm1, eax, 3

    // dest.argb = dest.argb + source.argb - dest.argb * sourec.alpha / 255
    paddw     mm0, mm1
    pmullw    mm1, mm2
    pmulhuw   mm1, mm6
    paddusw   mm1, mm5
    psrlw     mm1, 4
    psubsw    mm0, mm1
    pextrw    eax, mm0, 3

    // dest.rgb = dest.rgb * 255 / dest.alpha
    movq      mm1, mm0
    psllw     mm0, 8
    psubw     mm0, mm1
    pmulhuw   mm0, qword ptr MMDivTab[eax*8]
    packuswb  mm0, mm7
    movd      [edi], mm0
    mov       [edi].TARGBQuad.Alpha, al
end;

procedure CopyInterpolateData(var Deat: TImageData; const Source: TImageData; Alpha: Integer);
asm
    push      esi
    push      edi
    push      ebx
    cmp       ecx, 256
    jne       @@CvtPArgb
    cmp       [edx].TImageData.AlphaFlag, True
    je        @@CvtPArgb
    call      _SetCopyRegs
@@CpyLoop:
    push      ecx
    rep       movsd
    pop       ecx
    add       esi, eax
    add       edi, ebx
    dec       edx
    jne       @@CpyLoop
    jmp       @@Exit

@@CvtPArgb:
    push      ebp
    push      ecx
    call      _SetCopyRegs
    mov       ebp, eax
    call      SetMixerMM
    pop       eax
    movd      mm2, eax
    pshufw    mm2, mm2, 0
@@yLoop:
    push      ecx
@@xLoop:
    movd      mm0, [esi]
    punpcklbw mm0, mm7
    pshufw    mm1, mm0, 255
    pmullw    mm1, mm2
    psrlw     mm1, 8      // alpha0 = Source.Alpha * Alpha / 256
    pmullw    mm0, mm1    // Source.RGB = (Source.RGB * alpha0 + 127) / 255
    pmulhuw   mm0, mm6
    paddusw   mm0, mm5
    psrlw     mm0, 4
    packuswb  mm0, mm7
    movd      eax, mm1
    movd      [edi], mm0
    mov       [edi].TARGBQuad.Alpha, al
    add       esi, 4
    add       edi, 4
    loop      @@xLoop
    add       esi, ebp
    add       edi, ebx
    pop       ecx
    dec       edx
    jnz       @@yLoop
    pop       ebp
    emms
@@Exit:
    pop       ebx
    pop       edi
    pop       esi
end;

procedure FillBorder(var Data: TImageData; Radius: Integer; FillFlag: Integer);
asm
    push    ebp
    push    esi
    push    edi
    push    ebx
    test    ecx, 0fffh
    jne     @@yFill
    push    eax
    push    ecx
    mov     edi, [eax].TImageData.Stride
    imul    edi, edx
    add		  edi, [eax].TImageData.Scan0
    mov     ebp, [eax].TImageData.Width
    mov     ebx, [eax].TImageData.Height
    mov     esi, edx
    shl     esi, 1
    sub     ebx, esi
    sub     ebp, esi
    shl     ebp, 2
    shl     esi, 1
@@cLoop:
    mov     eax, [edi+esi]
    mov		  ecx, edx
    rep     stosd
    add     edi, ebp
    mov     eax, [edi-4]
    mov     ecx, edx
    rep		  stosd
    dec     ebx
    jnz     @@cLoop
    pop     ecx
    pop     eax
@@yFill:
    test    ecx, 0fff0000h
    jne     @@Exit
    mov     ebp, [eax].TImageData.Width
    mov		  edi, [eax].TImageData.Scan0
    mov     esi, [eax].TImageData.Stride
    imul    esi, edx
    add     esi, edi
    push    edx
@@tLoop:
    push    esi
    mov     ecx, ebp
    rep     movsd
    pop     esi
    dec     edx
    jnz     @@tLoop
    pop     edx
    mov     ebx, [eax].TImageData.Height
    sub     ebx, edx
    sub     ebx, edx
    imul    ebx, [eax].TImageData.Stride
    add     edi, ebx
    mov     esi, edi
    sub     esi, [eax].TImageData.Stride
@@bLoop:
    mov     ecx, ebp
    rep     movsd
    dec     edx
    jnz     @@bLoop
@@Exit:
    pop     ebx
    pop     edi
    pop     esi
    pop     ebp
end;

type
  TInterpolateProc = procedure;
  TElementsI = array[0..5] of Integer;

procedure _DoTransform(var Dest: TImageData; const Soutce: TImageData;
  const Elements: TElementsI; Radius: Integer; ipProc: TInterpolateProc);
var
  im22, im21, im12, im11: Integer;
  up, xDown, yDown: Integer;
  width, height, dstOffset: Integer;
  scan0: Pointer;
asm
    push    esi
    push    edi
    push    ebx
    mov     ebx, [ecx]
    mov     im11, ebx
    mov     ebx, [ecx+4]
    mov     im12, ebx
    mov     ebx, [ecx+8]
    mov     im21, ebx
    mov     ebx, [ecx+12]
    mov     im22, ebx
    push    dword ptr[ecx+16]
    push    dword ptr[ecx+20]
    mov     esi, Radius
    mov     ebx, [edx].TImageData.Width
    mov     edi, [edx].TImageData.Height
    sub     ebx, esi
    sub     edi, esi
    shl     ebx, 12
    shl     edi, 12
    mov     xDown, ebx          // xDown = (Source.Width - radius) * 4096
    mov     yDown, edi          // yDown = (Source.Height - radius) * 4096
    mov     ebx, [eax].TImageData.Height
    mov     height, ebx         // height = dest.Height
    mov     ebx, [eax].TImageData.Width
    mov     width, ebx          // width = dest.Width
    shl     ebx, 2
    neg     ebx
    add     ebx, [eax].TImageData.Stride
    mov     dstOffset, ebx      // dstOffset = dest.Stride - dest.Width * 4
    mov     edi, [eax].TImageData.Scan0
    mov     ebx, [edx].TImageData.Scan0
    mov     scan0, ebx
    mov     ebx, [edx].TImageData.Stride
    pop     ecx                 // ys
    pop     edx                 // xs
    shl     esi, 11
    mov     up, esi             // up = radius * 2048
    add     esi, 800h
    add     edx, esi            // xs += (up + 2048)
    add     ecx, esi            // ys += (up + 2048)
    pxor    xmm7, xmm7
    mov     eax, 40004h
    movd    xmm6, eax
    pshufd  xmm6, xmm6, 0
    call    SetMixerMM
@@yLoop:
    push    ecx                 // for (pd = (ARGB*)Dest.Scan0; height > 0; height --)
    push    edx                 // {
    push    width               //   x = xs
@@xLoop:                        //   y = ys
    cmp     ecx, up             //   for (w = width; w > 0; w --)
    jl      @@Next              //   {
    cmp     ecx, yDown
    jge     @@Next              
    cmp     edx, up
    jl      @@Next
    cmp     edx, xDown          //     if (y >= up && y < yDown && x >= up && x < xDown)
    jge     @@Next              //     {
    mov     esi, ecx            //       y0 = y / 4096
    mov     eax, edx            //       x0 = x / 4096
    sar     esi, 12
    sar     eax, 12
    imul    esi, ebx
    shl     eax, 2
    add     esi, eax
    add     esi, scan0          //       esi = Source.Scan0 + x0 * 4 + y0 * Source.Stride
    call    ipProc

    movd    eax, xmm0           //       xmm0 = eax = ipProc(src, x, y, esi)
    shr     eax, 24
    jz      @@Next
    call    MixerColor          //       if (eax >> 24) MixerColor(pd, mm0, eax)
@@Next:
    add     edi, 4              //       pd ++
    add     edx, im11           //       x += im11
    add     ecx, im12           //       y += im12
    dec     width               //     }
    jnz     @@xLoop
    pop     width
    pop     edx
    pop     ecx
    add     edi, dstOffset      //   (LPBYTE)pd += dstOffset
    add     edx, im21           //   xs += im21
    add     ecx, im22           //   ys += im22
    dec     height              // }
    jnz     @@yLoop
    emms
@@Exit:
    pop     ebx
    pop     edi
    pop     esi
end;

// 获取插值过程和扩展半径
function GetInterpolateProc(const Data: TImageData; var Proc: TInterpolateProc): Integer;
const
  procs: array[TInterpolateMode] of TInterpolateProc =
    (_GetBilinearColor, _GetNearColor, _GetBilinearColor, _GetBicubicColor);
  radius: array[TInterpolateMode] of Integer = (2, 1, 2, 4);
begin
  Proc := procs[Data.IpMode];
  Result := radius[Data.IpMode];
end;

// 按Matrix计算目标和源图像几何变换的裁剪矩形
function _GetTransformParams(dstWidth, dstHeight, srcWidth, srcHeight: Integer;
  var Matrix: TTransformMatrix; var dst, src: TRect): Boolean;

  function CalcRectI(Width, height: Integer; var r: TRect): Boolean;
  begin
    Inc(r.Right, r.Left);
    Inc(r.Bottom, r.Top);
    if r.Right > Width then r.Right := Width;
    if r.Bottom > Height then r.Bottom := Height;
    if r.Left > 0 then Dec(r.Right, r.Left) else r.Left := 0;
    if r.Top > 0 then Dec(r.Bottom, r.Top) else r.Top := 0;
    Result := (r.Right > 0) and (r.Bottom > 0);
  end;

var
  fx, fy, fwidth, fheight: Single;
begin
  Matrix.GetTransformSize(srcWidth, srcHeight, fx, fy, fwidth, fheight);
  Matrix.Invert;
  dst := Rect(Trunc(fx), Trunc(fy), _Infinity(fwidth), _Infinity(fheight));
  Result := CalcRectI(dstWidth, dstHeight, dst);
  if not Result then Exit;
  if (fx > 0) or (fy > 0) then
  begin
    if fx < 0 then fx := 0
    else if (fy < 0) then fy := 0;
    Matrix.Translate(fx, fy);
  end;
  Matrix.GetTransformSize(dst.Right, dst.Bottom, fx, fy, fwidth, fheight);
  src := Rect(Trunc(fx), Trunc(fy), _Infinity(fwidth), _Infinity(fheight));
  Result := CalcRectI(srcWidth, srcHeight, src);
  if not Result then Exit;
  if fx > 0 then Matrix.OffsetX := Matrix.OffsetX - fx;
  if fy > 0 then Matrix.OffsetY := Matrix.OffsetY - fy;
end;

procedure ImageTransform(var Dest: TImageData; x, y: Integer;
  const Source: TImageData; const Matrix: TTransformMatrix; Alpha: Single = 1);
var
  m: TTransformMatrix;
  e: TMatrixElements;
  eI: TElementsI;
  dstR, srcR: TRect;
  i, alphaI, radius: Integer;
  proc: TInterpolateProc;
  dst, src, tmp, sub: TImageData;
begin
  alphaI := Round(Alpha * 256);
  if alphaI <= 0 then Exit;
  if alphaI > 256 then alphaI := 256;
  m := TTransformMatrix.Create(Matrix);
  try
    m.OffsetX := m.OffsetX + x;
    m.OffsetY := m.OffsetY + y;
    if not _GetTransformParams(Dest.Width, Dest.Height,
      Source.Width, Source.Height, m, dstR, srcR) then  Exit;
    e := m.Elements;
    for i := 0 to 5 do
      eI[i] := Round(e.Elements[i] * 4096);
    radius := GetInterpolateProc(Source, proc);
    dst := GetSubImageData(Dest, dstR.Left, dstR.Top, dstR.Right, dstR.Bottom);
    tmp := GetSubImageData(Source, srcR.Left, srcR.Top, srcR.Right, srcR.Bottom);
    src := NewImageData(tmp.Width + radius * 2, tmp.Height + radius * 2);
    try
      src.AlphaFlag := Source.AlphaFlag;
      sub := GetSubImageData(src, radius, radius, tmp.Width, tmp.Height);
      CopyInterpolateData(sub, tmp, alphaI);
      FillBorder(src, radius, (eI[1] shl 16) or (eI[2] and $ffff));
      _DoTransform(dst, src, eI, radius, proc);
      if (eI[1] = 0) and (eI[2] = 0) and (alphaI = 256) and not Source.AlphaFlag
        and (Dest.Width = dst.Width) and (Dest.Height = dst.Height) then
        Dest.AlphaFlag := False;
    finally
      FreeImageData(src);
    end;
  finally
    m.Free;
  end;
end;

    代码很长，但核心代码还是_DoTransform过程和《Delphi图像处理 -- 平面几何变换类》中的三个插值过程。几何变换过程原理亦可参见该文。

    有关插值边界的处理，一般有2种办法，一是在插值过程中进行判断坐标是否超界而作相应的处理，二是舍弃边界部分，对于后者我是不主张的，因为那样是不完整的处理。我采用了扩展边框的办法进行边框插值处理，这样一来，虽然多了一道拷贝过程，却少了具体插值过程的坐标判断，二者抵消，插值速度应该是差不多的（据我测试，扩展边框办法在图像放大和旋转处理中速度还是略快一些），但是简化了插值代码。

    下面是图像几何变换的例子： 

procedure TForm1.Button1Click(Sender: TObject);
var
  bmp, newBmp: TBitmap;
  jpg: TJPEGImage;
  matrix: TTransformMatrix;
  source, dest: TImageData;
  r: TRect;
begin
  bmp := TBitmap.Create;
  matrix := TTransformMatrix.Create;
  try
    jpg := TJPEGImage.Create;
    try
      jpg.LoadFromFile('..\..\media\IMG_9440_mf.jpg');
      bmp.Assign(jpg);
    finally
      jpg.Free;
    end;
    matrix.Scale(1.2, 1.2);
    matrix.Shear(0.3, 0.3);
    r := matrix.GetTransformRect(bmp.Width, bmp.Height);
    if (r.Right <= 0) or (r.Bottom <= 0) then Exit;
    source := GetBitmapData(bmp);
    newBmp := TBitmap.Create;
    try
      newBmp.PixelFormat := pf32bit;
      newBmp.Width := r.Right;
      newBmp.Height := r.Bottom;
      dest := GetBitmapData(newBmp);
      SetInterpolateMode(source, imBicubic);
      ImageTransform(dest, 0, 0, source, matrix);
      Canvas.Draw(0, 0, newBmp);
    finally
      newBmp.Free;
    end;
  finally
    matrix.Free;
    bmp.Free;
  end;
end;

    代码中的SetInterpolateMode(source, imBicubic);语句为设置双立方插值方式，该函数也在《Delphi图像处理 -- 数据类型及公用过程》的ImageData.pas单元。

    运行效果截图如下：

    运行效果还是很好的，特别是边界效果。

    对于含Alpha信息的图像，必须作自乘预处理，这在CopyAlphaData函数中是自动判断后转换的。比较下面同一张png图片2种方式的线性插值几何变换处理效果截图，左边是使用自乘预处理后的插值处理效果，右边是没做自乘预处理后的插值处理效果：

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

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

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