Delphi图像处理 -- 图像卷积及高斯模糊

在图像的处理过程中，经常要用到卷积模板，如图像锐化、图像平滑、高斯模糊、Hough变换等，为此，本人使用Delphi编写了图像通用卷积处理过程和高斯模糊过程，代码如下：

过程定义：

  // 卷积图像。参数：
  //   Dest输出图，Source原图，Data自身操作图像
  //   ConvolTemplate卷积模板，必须为奇数矩阵
  //   Callback回调函数，返回True终止操作，CallbackData回调函数参数地址
  procedure ImageConvolution(var Data: TImageData; ConvolTemplate: array of Integer); overload;
  procedure ImageConvolution(var Dest: TImageData; const Source: TImageData;
    ConvolTemplate: array of Integer); overload;
  // 无效参数或者被回调函数终止操作返回False。
  function ImageConvolution(var Dest: TImageData; const Source: TImageData;
    ConvolTemplate: array of Integer;
    Callback: TImageAbort; CallbackData: Pointer): Boolean; overload;

实现代码：

procedure ConvolutionMMX(Template: Pointer; Size, TemplateOffset: Integer); pascal;
var
  width, height: Integer;
  dstOffset, srcOffset: Integer;
asm
    mov         width, ecx
    mov         height, edx
    mov         dstOffset, ebx
    mov         srcOffset, eax
    mov         eax, TemplateOffset
    pxor        mm7, mm7
@@yLoop:                        // for (y = 0; y < Dst.Height; y ++)
    push        width           // {
@@xLoop:                        //   for (x = 0; x < Dst.Width; x ++)
    push        esi             //   {
    pxor        mm0, mm0        //     mm0 = 0
    mov         ebx, Template
    mov         edx, Size       //     for (I = 0; I < TemplateSize; I ++)
@@iLoop:                        //     {
    mov         ecx, Size       //       for (J = 0; J <= TemplateSize; J ++)
@@jLoop:                        //       {
    movd        mm1, [esi]
    punpcklbw   mm1, mm7
    pmullw      mm1, [ebx]      //         mm1 = pixels[I, J] * Template[i * j];
    paddsw      mm0, mm1        //         mm0 += mm1
    add         esi, 4          //         esi += 4
    add         ebx, 8          //         edx += 8
    loop        @@jLoop         //       }
    add         esi, eax        //       esi += TemplateOffset
    dec         edx
    jnz         @@iLoop         //     }
    packuswb    mm0, mm7
    mov         dl, [edi].TARGBQuad.Alpha
    movd        [edi], mm0      //     *(ARGB*)edi = mm0
    mov         [edi].TARGBQuad.Alpha, dl
    pop         esi
    add         esi, 4          //     esi += 4
    add         edi, 4          //     edi += 4
    dec         width
    jnz         @@xLoop         //   }
    add         esi, srcOffset  //   esi += srcOffset
    add         edi, dstOffset  //   edi += dstOffset
    pop         width
    dec         height
    jnz         @@yLoop         // }
    emms
end;

procedure Convolution(Template: Pointer; Size, TemplateOffset: Integer); pascal;
var
  width, height: Integer;
  dstOffset, srcOffset: Integer;
  reds, greens, blues: Integer;
asm
    mov     width, ecx
    mov     height, edx
    mov     dstOffset, ebx
    mov     srcOffset, eax
@@yLoop:
    push    width
@@xLoop:
    push    esi
    push    edi
    mov     edi, Template
    xor     ebx, ebx
    mov     reds, ebx
    mov     greens, ebx
    mov     blues, ebx
    mov     ecx, Size
@@iLoop:
    push    ecx
    mov     ecx, Size
@@jLoop:
    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     esi, 4
    add     edi, 4
    loop    @@jLoop
    add     esi, TemplateOffset
    pop     ecx
    loop    @@iLoop
    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
    dec     width
    jnz     @@xLoop
    add     esi, srcOffset
    add     edi, dstOffset
    pop     width
    dec     height
    jnz     @@yLoop
end;

function ImageConvolution(var Dest: TImageData; const Source: TImageData;
  ConvolTemplate: array of Integer;
  Callback: TImageAbort; CallbackData: Pointer): Boolean;
var
  i, j, n: Integer;
  p: PInteger;
  Dst, Src: TImageData;
  Template: Pointer;
  TemplateSize, Radius: Integer;
  TemplateOffset: Integer;
  Nuclear: Integer;
  Proc: Pointer;
begin
  Result := False;
  n := Length(ConvolTemplate);
  if ImageEmpty(Dest) or ImageEmpty(Source) or (n = 0) then
    Exit;
  TemplateSize := Trunc(Sqrt(n));
  Radius := TemplateSize shr 1;
  GetExpandData(Dest, Source, Radius, Dst, Src);
  TemplateOffset := Src.Stride - (TemplateSize shl 2);
  Nuclear := 0;
  for i := 0 to n - 1 do
    Inc(Nuclear, ConvolTemplate[i]);
  if Nuclear <= 1 then
  begin
    GetMem(Template, n * Sizeof(int64));
    j := 0;
    for i := 0 to n - 1 do
    begin
      TWordArray(Template^)[j] := ConvolTemplate[i];
      TWordArray(Template^)[j + 1] := ConvolTemplate[i];
      TWordArray(Template^)[j + 2] := ConvolTemplate[i];
      TWordArray(Template^)[j + 3] := ConvolTemplate[i];
      Inc(j, 4);
    end;
    Proc := @ConvolutionMMX;
    if @Dest <> @Source then
      CopyAlpha(Dst, Source);
  end
  else
  begin
    GetMem(Template, n * Sizeof(Integer));
    p := Template;
    for i := 0 to n - 1 do
    begin
      p^ := ConvolTemplate[i] * 65536 div Nuclear;
      Inc(p);
    end;
    Proc := @Convolution;
  end;
  try
    if Assigned(Callback) then
      Result := ExecuteAbort(Dst, Src, Proc,
        [Template, TemplateSize, TemplateOffset], Callback, CallbackData)
    else
      Result := ExecuteProc(Dst, Src, Proc,
        [Template, TemplateSize, TemplateOffset]);
  finally
    FreeImageData(Src);
    FreeMem(Template);
  end;
end;

procedure ImageConvolution(var Data: TImageData; ConvolTemplate: array of Integer);
begin
  ImageConvolution(Data, Data, ConvolTemplate, nil, nil);
end;

procedure ImageConvolution(var Dest: TImageData; const Source: TImageData;
  ConvolTemplate: array of Integer);
begin
  ImageConvolution(Dest, Source, ConvolTemplate, nil, nil);
end;

因为要进行图像卷积操作，必须拷贝出源图像，本文过程在拷贝源图像时使用ImageGetExpandData过程附带扩展了其边框，这样在卷积图像边界像素时不必再进行坐标范围判断，从而简化了代码和提高了运行速度。

由于本文提供的是通用的卷积过程，尽管使用了BASM代码，但还是比较耗时的，如果对速度有特别要求的卷积操作，建议写专用的处理过程为好，例如Hough变换，其卷积核为1，就没必要在卷积过程中使用耗时的除法运算了。

下面是高斯模糊过程测试代码：

var
  jpg: TJPEGImage;
  Data: TImageData;
begin
  jpg := TJPEGImage.Create;
  jpg.LoadFromFile('D:\VclLib\GdiplusDemo\Media\20041001.jpg');
  Canvas.Draw(0, 0, jpg);
  Data := GetImageData(jpg);
  ImageGaussiabBlur(Data, 3);
  ImageDataAssignTo(Data, jpg);
  Canvas.Draw(0, jpg.Height, jpg);
  FreeImageData(Data);
  jpg.Free;
end;

下面是采用不同Q值进行高斯模糊的运行结果图：

左边是原始图像，其后分别是采用Q=1、2、3调用GdipGaussiabBlur过程高斯模糊后的图像，最右边一幅是使用Photoshop8.1cs高斯模糊（Q=3）处理的。

可以看出，本过程处理的图片和Photoshop处理的图片效果基本接近，同等Q值（3）情况下，Photoshop处理的模糊程度似乎稍稍强一些，可能是算法上的差异，我的过程自动计算的Radius是Q+2，我用Q=3，Radius=6，也就是采用13*13卷积矩阵的模糊效果似乎更接近Photoshop的处理效果，难道Photoshop的Radius为Q * 2？这样一来，采用的卷积矩阵会更大，应该不是，不过，上面的高斯模糊过程中Radius可以任意给定而不通过自动计算，可以反复试试，以模仿Photoshop的处理；另外，当Q=3时图像好像已经模糊的很厉害了，其实，对于高分辨率（200-300）图像，Q=3的高斯模糊程度并没有这么明显，我采用的图片分辨率为72，而且这张图片已经是经过Photoshop过滤、模糊处理后合成的。

下面给出一个直接调用GdipConvolution过程进行Hough变换的例子作为本文结束：

const
//  Ray: array[0..8] of Integer = (-1, 0, 1, -1, 0, 1, -1, 0, 1);
//  Ray: array[0..8] of Integer = (-1, -1, 0, -1, 0, 1, 0, 1, 1);
    Ray: array[0..8] of Integer = (-1, -1, -1, 0, 0, 0, 1, 1, 1);
var
  bmp: TGpBitmap;
  Data: TImageData;
  g: TGpGraphics;
begin
  bmp := TGpBitmap.Create('D:\msn.jpg');
  g := TGpGraphics.Create(Canvas.Handle);
  g.DrawImage(bmp, 0, 0);
  Data := GetImageData(bmp);
  ImageConvolution(Data, Ray);
  ImageDataAssignTo(Data, bmp, False);
  g.DrawImage(bmp, 0, Data.Height);
  g.Free;
  FreeImageData(Data);
  bmp.Free;
end;

运行效果如下图:

左边是原始图像，后面分别为用不同的Hough矩阵的运行效果，见例子代码中的几组Ray数组。作Hough变换，一般应该采用灰度图，本例子只是演示一下，没做图片灰度处理。

文章中所用数据类型及一些过程见《Delphi图像处理 -- 数据类型及内部过程》和《Delphi图像处理 -- 图像像素结构与图像数据转换》。

例子使用GDI+版本下载地址和说明见《GDI+ for VCL基础 -- GDI+ 与 VCL》。

尽管我十分努力，但水平有限，错误在所难免，欢迎指正和指导。邮箱地址：

[e-mail]

注：本文已经于2009.10.28重新整理。

本文代码于2010.5.20重新修订过。增加了拷贝形式的调整过程和响应回调函数的调整过程。代码中的ExecuteAbort过程和ExecuteProc过程见《Delphi图像处理 -- 图像像素结构与图像数据转换》。

修订过的卷积过程对于卷积核小于等于1的卷积模板采用了MMX处理，过程速度相当快，当然也要求其模板单个数据项不得超过127，整个模板数据项与255的乘积和不得超过$7fff（一般卷积核小于等于1的卷积模板都远远小于这些值）。对于卷积核大于1的卷积模板，本文代码也作了很大改进，速度比以前有所提高。

另外，高斯模糊处理过程已经另写（见《Delphi图像处理 -- 图像高斯模糊处理（改进版）》），本文的高斯模糊例子就作为遗迹存放吧（代码是有效的，只是处理过程不在本文）！