2013-08-18 20:18:22 GoodShot 阅读数 5666

图像处理中的卷积与模板

1.使用模板处理图像相关概念:      

     模板:矩阵方块,其数学含义是一种卷积运算。

  卷积运算:可看作是加权求和的过程,使用到的图像区域中的每个像素分别与卷积核(权矩阵)的每个元素对应相乘,所有乘积之和作为区域中心像素的新值。

     卷积核:卷积时使用到的权,用一个矩阵表示,该矩阵与使用的图像区域大小相同,其行、列都是奇数,是一个权矩阵。

     卷积示例:

             假设3 * 3的像素区域R与卷积核G分别为:


则卷积运算为:

              R5(中心像素)=R1G1 + R2G2 + R3G3 + R4G4 + R5G5 + R6G6 + R7G7 + R8G8 + R9G9

           

 

2、使用模板处理图像时涉及到的问题:

      边界问题:当处理图像边界像素时,卷积核与图像使用区域不能匹配,卷积核的中心与边界像素点对应,卷积运算将出现问题。

      处理办法:

              A.忽略边界像素,即处理后的图像将丢掉这些像素。

              B.保留原边界像素,即copy边界像素到处理后的图像。

 

3、常用模板:

      (a)低通滤波器

                             
 

   (b)高通滤波器                 

                       

   (c)平移和差分边缘检测

                      

    (d)匹配滤波边缘检测

                       

 

    (e)边缘检测

                           

 

    (f)梯度方向边缘检测

                      

 

                     

        

 

4、我们来看一下一维卷积的概念.

    卷积(convolution,另一个通用名称是德文的Faltung)的名称由来,是在于当初定义它时,定义成integ(f1(v)*f2(t-v))dv,积分区间在0到t之间。举个简单的例子,大家可以看到,为什么叫“卷积”了。比方说在(0,100)间积分,用简单的辛普生积分公式,积分区间分成100等分,那么看到的是f1(0)和f2(100)相乘,f1(1)和f2(99)相乘,f1(2)和f2(98)相乘,.........等等等等,就象是在坐标轴上回卷一样。所以人们就叫它“回卷积分”,或者“卷积”了。

    连续空间的卷积定义是f(x)与g(x)的卷积是f(t-x)g(x)在t从负无穷到正无穷的积分值.t-x要在f(x)定义域内,所以看上去很大的积分实际上还是在一定范围的.
       实际的过程就是f(x)先做一个Y轴的反转,然后再沿X轴平移t就是f(t-x),然后再把g(x)拿来,两者乘积的值再积分.想象一下如果g(x)或者f(x)是个单位的阶越函数.那么就是f(t-x)与g(x)相交部分的面积.这就是卷积了.

    卷积运算满足交换律,也就是说:f与g进行卷积完全等于g与f进行卷积。

   由两个函数f和g进行卷积而得到的函数f*g,一般要比原来的f和g都要光滑。所以在图像处理中对图像进行卷积后会使原图像模糊。因为卷积具有平滑作用。

    有趣的是,如果把两个人的照片互相进行卷积,所得到的照片,就同时和这两个人都很相像。

把积分符号换成求和就是离散空间的卷积定义了.

   那么在图像中卷积是什么意思呢,就是图像就是图像f(x),模板是g(x),然后将模版g(x)在模版中移动,每到一个位置,就把f(x)与g(x)的定义域相交的元素进行乘积并且求和,得出新的图像一点,就是被卷积后的图像.模版又称为卷积核.卷积核做一个矩阵的形状。

5、卷积运算时的核函数

      在Matlab中,对图像进行卷积运算时,都要先得到一个核函数,其实就是模板。其函数调用是:

>> G=fspecial('gaussian',5,0.5)

G =

     0.0000    0.0000    0.0002    0.0000    0.0000

    0.0000    0.0113    0.0837    0.0113    0.0000

    0.0002    0.0837    0.6187    0.0837    0.0002

    0.0000    0.0113    0.0837    0.0113    0.0000

     0.0000    0.0000    0.0002    0.0000    0.0000

>> G=fspecial('gaussian',5,1.5)

G =

    0.0144    0.0281    0.0351    0.0281    0.0144

   0.0281    0.0547    0.0683    0.0547    0.0281

   0.0351    0.0683    0.0853    0.0683    0.0351

   0.0281    0.0547    0.0683    0.0547    0.0281

   0.0144    0.0281    0.0351    0.0281    0.0144

    能够看出来,fspesial()函数的第一个参数表示返回高斯核函数(低通滤波器、模板等名称其实都一样)。第二个参数“5”表示该模板的大小,是5X5的矩阵。第三个参数是sigma了。

2018-06-07 15:23:46 zuange2417 阅读数 754

根据冈萨雷斯一书中的图像卷积概念。

卷积,是把核矩阵顺时针旋转180度,在和图像做变换。这个变换类似书中有,这里不过多解释。

我个人认为,虽然国外都是convolution,译为数学中的卷积是没错,并且人家就是这个数学上的卷积的概念

但其实根据卷积定义,用在图像——二维离散上,积分变为加法,其实更像是实现了矩阵的点积。

所以,当大部分情况,核矩阵旋转后不改变的情况下,当成点积更好理解。

2007-08-22 11:21:00 maozefa 阅读数 9126

阅读提示:

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

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

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

   本文代码必须包括文章《Delphi图像处理 -- 数据类型及公用过程》中的ImageData.pas单元

 

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

procedure Convolution32(var Dest: TImageData; const Source: TImageData;
  Template: Pointer; Size, TemplateOffset: Integer);
var
  width, height: Integer;
  dstOffset, srcOffset: Integer;
asm
    push      esi
    push      edi
    push      ebx
    push      ecx
    call      _SetCopyRegs
    mov       width, ecx
    mov       height, edx
    mov       dstOffset, ebx
    mov       srcOffset, eax
    mov       eax, TemplateOffset
    pop       ebx             // ebx = Template
    pxor      xmm7, xmm7
@@yLoop:                      // for (y = 0; y < Dst.Height; y ++)
    push      width           // {
@@xLoop:                      //   for (x = 0; x < Dst.Width; x ++)
    push      esi             //   {
    push      ebx
    pxor      xmm0, xmm0      //     xmm0 = 0
    mov       edx, size       //     for (I = 0; I < TemplateSize; I ++)
@@iLoop:                      //     {
    mov       ecx, size       //       for (J = 0; J <= TemplateSize; J ++)
@@jLoop:                      //       {
    movd      xmm1, [esi]
    punpcklbw xmm1, xmm7
    punpcklwd xmm1, xmm7
    cvtdq2ps  xmm1, xmm1
    mulps     xmm1, [ebx]     //         xmm1 = pixels[I, J] * Template[i * j]
    addps     xmm0, xmm1      //         xmm0 += xmm1
    add       esi, 4          //         esi += 4
    add       ebx, 16         //         edx += 16
    loop      @@jLoop         //       }
    add       esi, eax
    dec       edx
    jnz       @@iLoop         //     }
    cvtps2dq  xmm0, xmm0
    packssdw  xmm0, xmm7
    packuswb  xmm0, xmm7
    movd      [edi], xmm0     //     *(ARGB*)edi = xmm0
    pop       ebx
    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         // }
    pop       ebx
    pop       edi
    pop       esi
end;

procedure ConvolutionOne(var Dest: TImageData; const Source: TImageData;
  Template: Pointer; Size, TemplateOffset: Integer);
var
  width, height: Integer;
  dstOffset, srcOffset: Integer;
asm
    push      esi
    push      edi
    push      ebx
    push      ecx
    call      _SetCopyRegs
    mov       width, ecx
    mov       height, edx
    mov       dstOffset, ebx
    mov       srcOffset, eax
    mov       eax, TemplateOffset
    pop       ebx             // ebx = Template
    pxor      mm7, mm7
@@yLoop:                      // for (y = 0; y < Dst.Height; y ++)
    push      width           // {
@@xLoop:                      //   for (x = 0; x < Dst.Width; x ++)
    push      esi             //   {
    push      ebx
    pxor      mm0, mm0        //     mm0 = 0
    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
    dec       edx
    jnz       @@iLoop         //     }
    packuswb  mm0, mm7
    movd      [edi], mm0      //     *(ARGB*)edi = mm0
    pop       ebx
    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
    pop       ebx
    pop       edi
    pop       esi
end;

procedure CvtTemplate8(Template8, Template: Pointer; n: Integer);
asm
    push      eax
    push      ecx
    pcmpeqw   mm7, mm7
    psrlq     mm7, 16
@@Loop:
    movd      mm0, [edx]
    pshufw    mm0, mm0, 0
    pand      mm0, mm7
    movq      [eax], mm0
    add       edx, 4
    add       eax, 8
    loop      @@Loop
    pop       ecx
    pop       eax
    shr       ecx, 1
    shl       ecx, 3
    add       eax, ecx
    add       eax, 6
    mov       word ptr[eax], 1  // (short)Template8[n / 2][3] = 1
    emms
end;

procedure CvtTemplate16(Template16, Template: Pointer; n: Integer);
var
  nuclear: Single;
asm
    push      edx
    push      ecx
    fldz
@@SumLoop:
    fadd      dword ptr[edx]
    add       edx, 4
    loop      @@SumLoop
    fstp      nuclear
    fwait
    pop       ecx
    pop       edx
    cmp       nuclear, 0
    je        @@1
    movss     xmm7, nuclear
    shufps    xmm7, xmm7, 0
@@Loop:
    movss     xmm0, [edx]
    shufps    xmm0, xmm0, 0
    divps     xmm0, xmm7
    movaps    [eax], xmm0
    add       edx, 4
    add       eax, 16
    loop      @@Loop
    jmp       @@Exit
@@1:
    push      eax
    push      ecx
    pcmpeqd   xmm7, xmm7
    psrldq    xmm7, 4         // xmm7 >> 32
@@Loop1:
    movss     xmm0, [edx]
    shufps    xmm0, xmm0, 0
    andps     xmm0, xmm7
    movaps    [eax], xmm0
    add       edx, 4
    add       eax, 16
    loop      @@Loop1
    pop       ecx
    pop       eax
    shr       ecx, 1
    shl       ecx, 4
    add       eax, ecx
    add       eax, 12
    fld1
    fstp      dword ptr[eax]  // (float)Template16[n / 2][3] = 1.0
@@Exit:
end;

procedure CvtTemplate16I(Template16, Template: Pointer; n, nuclear: Integer);
asm
    cmp       nuclear, 1
    jle       @@1
    movd      xmm7, nuclear
    pshufd    xmm7, xmm7, 0
    cvtdq2ps  xmm7, xmm7
@@Loop:
    movd      xmm0, [edx]
    pshufd    xmm0, xmm0, 0
    cvtdq2ps  xmm0, xmm0
    divps     xmm0, xmm7
    movaps    [eax], xmm0
    add       edx, 4
    add       eax, 16
    loop      @@Loop
    jmp       @@Exit
@@1:
    push      eax
    push      ecx
    pcmpeqd   xmm7, xmm7
    psrldq    xmm7, 4         // xmm7 >> 32
@@Loop1:
    movd      xmm0, [edx]
    pshufd    xmm0, xmm0, 0
    pand      xmm0, xmm7
    cvtdq2ps  xmm0, xmm0
    movaps    [eax], xmm0
    add       edx, 4
    add       eax, 16
    loop      @@Loop1
    pop       ecx
    pop       eax
    shr       ecx, 1
    shl       ecx, 4
    add       eax, ecx
    add       eax, 12
    fld1
    fstp      dword ptr[eax]  // (float)Template16[n / 2][3] = 1.0
@@Exit:
end;

procedure ImageConvolutionI(var Data: TImageData; const Template: array of Integer);
var
  i, n, sum: Integer;
  buffer, pTemplate: Pointer;
  templateSize, templateOffset: Integer;
  nuclear: Integer;
  isFloatOp: Boolean;
  src: TImageData;
begin
  if Data.AlphaFlag then
    ArgbConvertPArgb(Data);
  n := Length(Template);
  templateSize := Trunc(Sqrt(n));
  src := _GetExpandData(Data, templateSize shr 1);
  templateOffset := src.Stride - (templateSize shl 2);
  nuclear := 0;
  for i := 0 to n - 1 do
    Inc(nuclear, Template[i]);
  isFloatOp := nuclear > 1;
  if not isFloatOp then
  begin
    sum := 0;
    for i := 0 to n - 1 do
    begin
      Inc(sum, Template[i] * 255);
      isFloatOp := (sum > 32767) or (sum < -32768);
      if isFloatOp then Break;
    end;
  end;
  if not isFloatOp then
  begin
    buffer := GlobalAllocPtr(GHND, n * Sizeof(TMMType));
    CvtTemplate8(buffer, @Template[0], n);
    ConvolutionOne(Data, src, buffer, templateSize, templateOffset);
  end
  else
  begin
    buffer := GlobalAllocPtr(GHND, n * 16 + 12);
    pTemplate := Pointer((Integer(buffer) + 12) and (not 15));
    CvtTemplate16I(pTemplate, @Template[0], n, nuclear);
    Convolution32(Data, src, PTemplate, templateSize, templateOffset);
  end;
  GlobalFreePtr(Buffer);
  FreeImageData(src);
  if Data.AlphaFlag then
    PArgbConvertArgb(Data);
end;

procedure ImageConvolution(var Data: TImageData; const Template: array of Single);
var
  n: Integer;
  buffer, pTemplate: Pointer;
  templateSize, templateOffset: Integer;
  src: TImageData;
begin
  if Data.AlphaFlag then
    ArgbConvertPArgb(Data);
  n := Length(Template);
  templateSize := Trunc(Sqrt(n));
  src := _GetExpandData(Data, templateSize shr 1);
  templateOffset := src.Stride - (templateSize shl 2);
  buffer := GlobalAllocPtr(GHND, n * 16 + 12);
  pTemplate := Pointer((Integer(buffer) + 12) and (not 15));
  CvtTemplate16(pTemplate, @Template[0], n);
  Convolution32(Data, src, PTemplate, templateSize, templateOffset);
  GlobalFreePtr(Buffer);
  FreeImageData(src);
  if Data.AlphaFlag then
    PArgbConvertArgb(Data);
end;

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

   考虑不同的应用环境,本文写了浮点卷积矩阵和整数卷积矩阵2个版本。

   浮点版本采用SSE浮点运算,每次可同时处理图像像素ARGB4个分量,运行速度比较快;

   整数版本有2种处理方法,对于卷积核小于等于1的卷积模板采用了MMX处理,过程速度相当快,当然也要求其模板单个数据项不得超过127,整个模板数据项与255的乘积和不得超过$7fff,一般卷积核小于等于1的卷积模板都远远小于这些值(当然,代码中已经对此作了检查,并做出了相应的处理)。对于卷积核大于1的卷积模板,则转换为浮点数版本处理。 

        下面给出一个直接调用GdipConvolution过程进行Hough变换的例子:

procedure TForm1.Button3Click(Sender: TObject);
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;
  g: TGpGraphics;
  data: TImageData;
begin
  bmp := TGpBitmap.Create('..\media\msn.jpg');
  g := TGpGraphics.Create(Canvas.Handle);
  g.DrawImage(bmp, 0, 0);
  data := LockGpBitmap(bmp);
  ImageConvolutionI(Data, Ray);
  UnlockGpBitmap(bmp, data);
  g.DrawImage(bmp, data.Width, 0);
  g.Free;
  bmp.Free;
end;

    运行效果如下图:

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

 

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

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

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

 

2017-07-28 13:50:23 TchChan 阅读数 1998

模板:矩阵方块,其数学含义是一种卷积运算

卷积运算:可看做是加权求和的过程,使用到的图像区域中的每个像素分别于卷积核(权矩阵)的每个元素对应相乘,所有乘机之和作为区域中心像素的新值。

利用卷积可以实现对图像模糊处理,边缘检测,产生轧花效果的图像。

卷积核:卷积时用到的权用一个矩阵表示,该矩阵与使用的图像区域大小相同,其行列都是奇数(这一点很重要)

举个例子:

一个3*3的像素区域R与卷积核G的卷积运算:

R5(中心像素)=R1G1+R2G2+R3G3+R4G4+R5G5+R6G6+R7G7+R8G8+R9G9


下篇将说明一下卷积运算遇到的问题和几种处理办法以及模板






2017-10-21 19:21:45 tsfx051435adsl 阅读数 3625

1、一维信号的卷积
  卷积是一种是数学运算,对于一维连续信号f(x)来说,f(x)和g(x)的卷积是一种积分运算,如下
  这里写图片描述
  它的几何意义是这样的:第一步,先把g(x)关于原点对称,得到g(-x);第二步,把g(-x)向右平移a个单位得到g(a-x)的图像,如果两个图像有交叉部分,那么交叉部分的面积就是卷积。
  上面说的是一维连续信号,如果是一维离散信号呢,因为积分是求和的极限,那么对于离散信号的卷积就是一种求和运算。对于离散信号x(n)和h(n)来说,他们的卷积如下
  这里写图片描述
  几何过程也是类似的。先关于远点对称,在平移,只不过只能平移整数个单位,因为是离散信号,再相乘再相加。
2、一维信号的卷积有什么用
  对于信号处理来说,基本的模型是输入信号,处理信号,输出信号这个过程,不管什么信号,总得有输出,从输入到输出的过程都叫信号处理过程,其中包含的噪声、滤波等等都是这个过程中发生的事情。可以证明(谁证明的就不知道了),假设输入信号叫x,处理信号对应的系统函数叫h,那么输入信号等于x和h的卷积。(具体的证明见信号和系统这本书)。
  但是卷积是积分运算,计算起来比较麻烦,所以能不能不通过卷积运算就能把输出信号算出来呢?能!就是通过变换域的方法,什么叫变换域,在正常处理中我们的信号都是关于时间t的信号,也就是说是自变量是t,我们叫时间域(就类似于定义域)。通过一种可逆的变换把时间域变换到别的域(比如傅里叶变换,拉普拉斯变换),从别的域进行处理,再反变换(因为我们需要的是时间域的东西,所以必须得可逆)回来。傅里叶变换和拉普拉斯变换分别是把时间域(简称时域)变换到频率域和s域(就是指变换之后的自变量是频率和s),可以证明(证明过程见信号和系统书吧):时域的卷积相当于频域(s域)的相乘。
  所以,对于信号的处理可以在时域处理,也可以在变换域处理。
3、数字图像的卷积
  数字图像在实质上就是一个二维数组,在信号的角度上来说是二维离散信号,只不过自变量不是时间,而是空间像素点了,所以正常需要处理的是空间域,类比二维信号也有变换域(频域,小波域等等),对于二维信号的处理也有如下结论:处理后的图像g(x,y)等于输入图像f(x,y)和处理过程h(x,y)的卷积。其中A矩阵和B矩阵的卷积如下。
  这里写图片描述
  数字图像进行卷积干吗用?主要用来图像增强和去噪。
4、几个常用术语
  模板:模板就是一个矩阵(其实就是下面说的卷积核),也就是处理图像这个过程对应的函数,对应卷积运算。
  卷积运算:可以看做加权求和的过程,为了图像增强和减少噪声,用该像素点的邻域像素点进行加权,得到处理完之后这个点的像素点(会把噪声减弱,当然也会有别的问题)。
  卷积核:上面不是说了卷积运算了嘛,这个点的邻域和哪个矩阵进行加权啊?换句话说,这个点的邻域那么多像素点,这些像素点所占的权重是多少啊?这些权重构成一个矩阵,这个矩阵就叫卷积核,卷积核的行数和列数都是奇数。
  比如要对R5这个像素点进行卷积运算,有下面的式子成立,这个卷积核的数组直接决定卷积的结果。
  这里写图片描述
  使用卷积核的时候回出现一个问题——边界问题:当处理图像边界像素时,卷积核与图像使用区域不能匹配,卷积核的中心与边界像素点对应, 卷积运算将出现问题。
  处理办法:
A. 忽略边界像素,即处理后的图像将丢掉这些像素。
B. 保留原边界像素,即copy边界像素到处理后的图像。

5、常用的模板
这里写图片描述
这里写图片描述

卷积功能的实现

阅读数 105

没有更多推荐了,返回首页