2016-11-07 16:12:12 LG1259156776 阅读数 8882

图像特征提取:图像的矩特征

1. 矩的概念

图像识别的一个核心问题是图像的特征提取,简单描述即为用一组简单的数据(图像描述量)来描述整个图像,这组数据越简单越有代表性越好。良好的特征不受光线、噪点、几何形变的干扰。图像识别发展几十年,不断有新的特征提出,而图像不变矩就是其中一个。

矩是概率与统计中的一个概念,是随机变量的一种数字特征。设X为随机变量,c为常数,k为正整数。则量E[(xc)k]称为X关于c点的k阶矩。

比较重要的有两种情况:

1. c=0。这时ak=E(Xk)称为Xk阶原点矩

2. c=E(X)。这时μk=E[(XEX)k]称为Xk阶中心矩。

一阶原点矩就是期望。一阶中心矩μ1=0,二阶中心矩μ2就是X的方差Var(X)。在统计学上,高于4阶的矩极少使用。μ3可以去衡量分布是否有偏。μ4可以去衡量分布(密度)在均值附近的陡峭程度如何。

针对于一幅图像,我们把像素的坐标看成是一个二维随机变量(X,Y),那么一幅灰度图像可以用二维灰度密度函数来表示,因此可以用矩来描述灰度图像的特征。

不变矩(Invariant Moments)是一处高度浓缩的图像特征,具有平移、灰度、尺度、旋转不变性。M.K.Hu在1961年首先提出了不变矩的概念。1979年M.R.Teague根据正交多项式理论提出了Zernike矩。下面主要介绍这两种矩特征的算法原理与实现。

2. Hu矩

一幅M×N的数字图像f(i,j),其p+q阶几何矩mpq和中心矩μpq为:

 

mpq=i=1Mj=1Nipjqf(i,j)

 

 

μpq=i=1Mj=1N(ii¯)p(jj¯)qf(i,j)

 

其中f(i,j)为图像在坐标点(i,j)处的灰度值。i¯=m10/m00,j¯=m01/m00

若将m00看作是图像的灰度质量,则(i¯,j¯)为图像的质心坐标,那么中心矩μpa反映的是图像灰度相对于其灰度质心的分布情况。可以用几何矩来表示中心矩,0~3阶中心矩与几何矩的关系如下:

μ00=Mi=1Nj=1(ii¯)0(jj¯)0f(i,j)=m00

μ10=Mi=1Nj=1(ii¯)1(jj¯)0f(i,j)=0

μ01=Mi=1Nj=1(ii¯)0(jj¯)1f(i,j)=0

μ11=Mi=1Nj=1(ii¯)1(jj¯)1f(i,j)=m11y¯m10

μ20=Mi=1Nj=1(ii¯)2(jj¯)0f(i,j)=m20y¯m01

μ02=Mi=1Nj=1(ii¯)0(jj¯)2f(i,j)=m02y¯m01

μ30=Mi=1Nj=1(ii¯)3(jj¯)0f(i,j)=m302x¯m20+2x¯2m10

μ12=Mi=1Nj=1(ii¯)1(jj¯)2f(i,j)=m122y¯m11x¯m02+2y¯2m10

μ21=Mi=1Nj=1(ii¯)2(jj¯)1f(i,j)=m212x¯m11y¯m20+2x¯2m01

μ03=Mi=1Nj=1(ii¯)0(jj¯)3f(i,j)=m032y¯m02+2y¯2m01

为了消除图像比例变化带来的影响,定义规格化中心矩如下:

 

ηpq=μpaμγ00,(γ=p+q2,p+q=2,3,)

 

利用二阶和三阶规格中心矩可以导出下面7个不变矩组(Φ1 Φ7),它们在图像平移、旋转和比例变化时保持不变。

Φ1=η20+η02

Φ2=(η20η02)2+4η211

Φ3=(η203η12)2+3(η21η03)2

Φ4=(η30+η12)2+(η21+η03)2

Φ5=(η30+3η12)(η30+η12)[(η30+η12)23(η21+η03)2]+(3η21η03)(η21+η03)[3(η30+η12)2(η21+η03)2]

Φ6=(η20η02)[(η30+η12)2(η21+η03)2]+4η11(η30+η12)(η21+η03)

Φ7=(3η21η03)(η30+η12)[(η30+η12)23(η21+η03)2]+]+(3η12η30)(η21+η03)[3(η30+η12)2(η21+η03)2]

3. 利用OpenCV计算Hu矩

opencv里对Hu矩的计算有直接的API,它分为了两个函数:moments()函数用于计算中心矩,HuMoments函数用于由中心矩计算Hu矩。

Moments moments(InputArray array, bool binaryImage=false )

参数说明

  • 输入参数:array是一幅单通道,8-bits的图像,或一个二维浮点数组(Point of Point2f)。binaryImage用来指示输出图像是否为一幅二值图像,如果是二值图像,则图像中所有非0像素看作为1进行计算。
  • 输出参数:moments是一个类:
class Moments
{
public:
    Moments();
    Moments(double m00, double m10, double m01, double m20, double m11,
            double m02, double m30, double m21, double m12, double m03 );
    Moments( const CvMoments& moments );
    operator CvMoments() const;

里面保存了图像的2阶与3阶中心矩的值。

void HuMoments(const Moments& moments, double* hu)

参数说明:

  • 输入参数:moments即为上面一个函数计算得到的moments类型。
  • 输出参数:hu是一个含有7个数的数组。

int main(int argc, char** argv) 

    Mat image = imread(argv[1]);  
    cvtColor(image, image, CV_BGR2GRAY); 
    Moments mts = moments(image); 
    double hu[7]; 
    HuMoments(mts, hu); 
    for (int i=0; i<7; i++) 
    { 
        cout << log(abs(hu[i])) <<endl; 
    } 
   return 0; 
}

上面代码中,最终输出的值为log|Φi|

我们分别计算一幅图像在,旋转,噪声与模糊时的Hu矩。

qi qi qi qi

类别 log|Φ1| log|Φ2| log|Φ3| log|Φ4| log|Φ5| log|Φ6| log|Φ7|
原图 -6.76181 -19.1286 -23.7441 -26.776 -51.7618 -35.8491 -51.534
旋转 -6.72102 -19.0844 -23.5756 -25.9122 -51.4619 -35.4595 -50.7674
加放噪点 -6.76086 -19.1255 -23.7611 -26.3228 -51.5056 -35.895 -51.6321
模糊 -6.76183 -19.1295 -23.7451 -26.2767 -51.765 -35.8484 -51.5307

4. Zernike矩

Hu矩在图像描述上有广泛的应用,但是其低阶几何矩与图像整体特征有关,不包含太多的图像细节信息,而高阶几何矩易受噪声影响,因此很难利用几何矩恢复图像。

Zernike矩能够很容易地构造图像的任意高阶矩,并能够使用较少的矩来重建图像。Zernike矩是基于Zernike多项式的正交化函数,虽然其计算比较复杂,但是Zernide矩在图像旋转和低噪声敏感度方面具有较大的优越性。由于Zernike矩具有图像旋转不变性,而且可以构造任意高阶矩,所以被广泛应用对目标进行识别中。

4.1 Zernike矩多项式

首先要弄清楚什么是正交多项式。若函数W(x)在区间(a,b)可积,且W(x)0,则可作为权函数。

对于一个多项式的序列fi和权函数W(x),定义内积:<fm,fn>=bafm(x)fn(x)W(x)dx

nm,<fm,fn>=0,这些多项式则称为正交多项式。若fi除了正交之外,更有<fm,fn>=1的话,则称为规范正交多项式。

那么正交多项式有什么作用呢?答案是:逼近!正交多项式相当于基,任何一个n维多项式函数f(x)都可以用一组正交多项式加权求和来逼近。

 

Zernike在1934年提出了在单位圆上定义的一组正交多项式,即Zernike正交多项式,其定义形式为:

 

Rnm(ρ)=s=0(n|m|)/2(1)s[(ns)!]ρn2ss!(n+|m|2s)!(n+|m|2+s)!

 

 

Vnm(x,y)=Vnm(ρ,θ)=Rnm(ρ)ejmθ

 

其中Rnm(ρ)表示点(x,y)的径向多项式,Vnm(x,y)为Zernike正交多项式,n,m为正交多项式的阶数,n是非负整数,n|m|是偶数,并且n|m|

Zernike多项式Vnm(x,y)=Vnm(ρ,θ)是定义在单位圆x2+y21上的正交复函数的集合,具有重要的递推性质,即Rnm可由R(n2)mR(n4)m得到,公式如下:

 

Rnm(ρ)=[(K22ρ2+K3)R(n2)m(ρ)+K4R(n4)m(ρ)]K1

 

 

Rmm(ρ)=ρm

 

式中:K1=(n+1)(n1)(n2)/2,K2=2n(n1)(n2),K3=(n1)3,K4=n(n1)(n3)/2

4.2 Zernike矩的定义

由于Zernike多项式的正交完备性,所以在单位圆内的任何图像f(x,y)都可以唯一的用下面式子展开:

 

f(x,y)=n=0m=0ZnmVn,m(ρ,θ)

 

上式中的Znm就是Zernike矩。

对二维函数f(x,y)的Zernike矩的定义如下:

 

Znm=n+1π102π0[Vnm(ρ,θ)]f(ρ,θ)ρdydxdρdθ

 

 

=n+1πRnm(ρ)ejmθf(ρ,θ)dρdθ

 

式中ρ=x2+y2−−−−−−√(1<x,y<1)θ为轴xρ矢量在逆时针方向的夹角;Rnm(ρ)表示点(x,y)的径向多项式。

4.3 Zernike矩的计算

从Zernike矩的计算公式上来看,对于二维图像,其Zernike矩Znm为复数,将其实部和虚部分别记为CnmSnm,则有:

 

Cnm=2n+2π102π0[Rnm(ρ)cos(mθ)f(ρ,θ)ρdρdθ

 

 

Cnm=2n+2π102π0[Rnm(ρ)sin(mθ)f(ρ,θ)ρdρdθ

 

因为数字图像是离散形式的点,所以需要将上式离散化,把积分号换为求和号,但是需要作一些坐标变换。

对于N×N的图像f(x,y),令坐标原点位于图像的中心,则N/2x,yN/2,对于像素(x,y),引入2个参数(r,σ),唯一对应于像素,其定义为:

r=max(|x|,|y|)

如果|x|=r,则:

 

σ=2(rx)y|y|+xyr

 

如果|y|=r,则:

 

σ=2yxyr

 

我们容易计算出,r的取值范围为1N/2σ的取值范围是18r,再根据参数(r,σ)定义相应的极坐标:

 

ρ=2r/N,θ=πσ(4r)

 

所以,最终我们得到离散化的Zernike矩的计算公式:

 

Cnm=2n+2N2r=1N/2Rnm(2r/N)σ=18rcosπmσ4rf(r,σ)

 

 

Snm=2n+2N2r=1N/2Rnm(2r/N)σ=18rsinπmσ4rf(r,σ)

 

1. 确定图像的大小N×N,即公式中的N

2. 确定r,σ的范围;

3. 利用Zernike多项式的递推性质计算各阶Rnm(ρ),并结合上面Zernike矩计算公式,算出Cnm,Snm

4. 对Cnm,Snm求模,进而计算得到|Znm|

现在我们用Zernike矩来计算美女图像在4种状态下的值:

类别 log|Z11| log|Z20| log|Z22| log|Z31| log|Z40| log|Z42| log|Z44|
原图 11.1732 13.8469 12.3515 12.4391 14.2782 12.6137 11.5745
旋转 12.3036 13.8309 13.5861 12.0467 13.1320 13.8396 12.7862
加放噪点 11.1538 13.8490 12.3315 12.4316 14.2730 12.5925 11.5591
模糊 11.1636 13.8465 12.3480 12.4367 14.2799 12.6130 11.5752

通过表中,可以看出,Zernike在总体上效果比Hu矩更好(PS:感觉在旋转上好像差强人意!)

下面是Zernike矩的matlab实现[来自《现代数字图像-处理技术提高及应用案例详解》],这里偷懒了,有机会的话会把C++版的实现补上。

 View Code

5. 总结

不变矩的应用过程一般包括:

  1. 选择合适的不变矩类型;
  2. 选择分类器(如神经网络、最短距离等);
  3. 如果是神经网络分类器,则需要计算学习样例的不变矩去训练神经网络;
  4. 计算待识别对象的不变矩,输入神经网络就可得到待识别对象的类型,或者计算待识别对象不变矩与类别对象不变矩之间的距离,选择最短距离的类别作为待识别对象的类别。

可以看出,不变矩作用主要目的是描述事物(图像)的特征。人眼识别图像的特征往往又表现为“求和”的形式,因此不变矩是对图像元素进行了积分操作。

不变矩能够描述图像整体特征就是因为它具有平移不变形、比例不变性和旋转不变性等性质。

然而,另一方面图像的各阶不变矩究竟代表的什么特征很难进行直观的物理解释。

6. 参考资料

[1] 《现代数字图像处理》(matlab版)

[2] 正交多项式WIKI

[3] opencv形态描述

2018-12-10 10:00:49 yql_617540298 阅读数 240

一、图像

图像:对应矩阵,每个位置的元素就是该处的像素;

图像矩:一个从数字图形中计算出来的矩集,通常描述了该图像全局特征,并提出了大量的关于该图像不同类型的几何特征信息;(如大小、位置、方向和形状等);

一阶矩:与形状有关;

二阶矩:显示了曲线围绕直线平均值的扩展程度;

三阶矩:关于平均值的对称性测量;

通过二阶矩计算图像中心:python

import imutils
import cv2
import os

#ap = argparse.ArgumentParser()
#--image 参数: 磁盘中待处理图像的路径
#ap.add_argument("-i","--image",required=True,help="path to the input image")
#args = vars(ap.parse_args())

#image = cv2.imread(args["image"])
filepath = 'C:/Users/Administrator/Desktop/result/'
resultpath = 'C:/Users/Administrator/Desktop/result_zx/'
count = os.listdir(filepath)
for i in range(0,len(count)):
    print("i=",i)
    image = cv2.imread(filepath+str(i)+'.png')
    gray = cv2.cvtColor(image,cv2.COLOR_BGR2GRAY)
    #5×5 内核的高斯平滑
    blurred = cv2.GaussianBlur(gray,(5,5),0)
    #阈值化
    thresh = cv2.threshold(blurred,60,255,cv2.THRESH_BINARY)[1]
    
    #使用轮廓检测去定位这些白色区域
    #返回图像上每一个白块对应的边界点集合(即轮廓)
    cnts = cv2.findContours(thresh.copy(),cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_SIMPLE)
    cnts = cnts[0] if imutils.is_cv2() else cnts[1]
    
    #处理每一条轮廓
    for c in cnts:
        #计算轮廓区域图像的矩
        M = cv2.moments(c)
        print("M:",M)
        if M["m00"]==0:
            print("skip!")
        else:
            #计算轮廓的中心
            cX = int(M["m10"] / M["m00"])
            cY = int(M["m01"] / M["m00"])
            #在形状的中心 (cX, cY) 处绘制一个白色的小圆
            cv2.drawContours(image, [c], -1, (0, 255, 0), 2)
            cv2.circle(image, (cX, cY), 7, (255, 255, 255), -1)
            cv2.putText(image, "", (cX - 20, cY - 20),
                        cv2.FONT_HERSHEY_SIMPLEX, 0.5, (255, 255, 255), 2)
    cv2.imwrite(resultpath+str(i)+'.png',image)
    #cv2.imshow("Image", image)
    cv2.waitKey(0)

通过二阶矩计算图像重心:MATLAB

for i=0:29
    %I = imread('picture\0.png'); 
    bgFile = ['picture\',int2str(i),'.png'];% 读入图片的完整路径
    I = imread(bgFile);
    I = rgb2gray(I); 
    imshow(I); 
    I = double(I); 
    [rows,cols] = size(I); 
    x = ones(rows,1)*[1:cols]; 
    y = [1:rows]'*ones(1,cols); 
    area = sum(sum(I)); 
    meanx = sum(sum(I.*x))/area; 
    meany = sum(sum(I.*y))/area; 
    hold on; 
    plot(meanx,meany,'r+'); %十字标出重心位置
    F=getframe(gcf); % 获取整个窗口内容的图像
    %imwrite(F.cdata,'result\1.png')
    img = ['result\',num2str(i),'.png']
    imwrite(F.cdata,img); 
    hold off;
end;

二、中心距求主轴方向

引用原文链接:https://blog.csdn.net/hong__fang/article/details/49851569

1. 中心距

2. 中心距求主轴方向

主轴方向与x轴正向的最小夹角,如果角度位于y正方向,则,位于y负方向,

求中心距求主轴另外公式:

求中心距时,如果是二值图像,可以只根据图像边界坐标求主方向,与利用图像所有坐标求主轴方向相比,只利用边界坐标会存在一定偏差,约

MATLAB实现:主轴方向

function test
 
    I = imread('picture\0.png');
  
    [cm ju] = qijieju(uint8(I));
    m00 = cm(1);
    mu11 = cm(2);
    mu02 = cm(3);
    mu20 = cm(4);
    a = mu20 / m00; 
    b = mu11 / m00;
    c = mu02 / m00;
    square = sqrt( 4 * b * b + (a - c) * (a - c) );
	%求主轴方法1
    theta1 = atan2( 2 * b, a - c + square )*180/pi
    %save('output\output.txt','theta1','-ascii')
    %B = imrotate(I,theta1,'bicubic')
    %imwrite(B,'output\1.png')
    %求主轴方法2
    theta2 = atan2(2*mu11,(mu20-mu02))/2*180/pi
    I2 = imrotate(I,-theta2*10,'nearest')
    imwrite(I2,'output\0.png')
end
 
%求不变矩及中心矩
function [cm ju] = qijieju(I0)
    A=double(I0);
    [nc,nr]=size(A);
    [x,y]=meshgrid(1:nr,1:nc);
    x=x(:);
    y=y(:);
    A=A(:);
    m00=sum(A);
    if m00==0
        m00=eps;
    end
    m10=sum(x.*A);
    m01=sum(y.*A);
    xmean=m10/m00; %重心
    ymean=m01/m00;
    cm00=m00; %归一化中心矩
    cm02=(sum((y-ymean).^2.*A))/(m00^2);
    cm03=(sum((y-ymean).^3.*A))/(m00^2.5);
    cm11=(sum((x-ymean).*(y-ymean).*A))/(m00^2);
    cm12=(sum((x-ymean).*(y-ymean).^2.*A))/(m00^2.5);
    cm20=(sum((x-xmean).^2.*A))/(m00^2);
    cm21=(sum((x-xmean).^2.*(y-ymean).*A))/(m00^2.5);
    cm30=(sum((x-xmean).^3.*A))/(m00^2.5);
    ju(1)=cm20+cm02;  %
    ju(2)=(cm20-cm02)^2+4*cm11^2; %
    ju(3)=(cm30-3*cm12)^2+(3*cm21-cm03)^2; %
    ju(4)=(cm30+cm12)^2+(cm21+cm03)^2;  %
    ju(5)=(cm30-3*cm12)*(cm30+cm12)*((cm30+cm12)^2-3*(cm21+cm03)^2)+(3*cm21-cm03)*(cm21+cm03)*(3*(cm30+cm12)^2-(cm21+cm03)^2); %
    ju(6)=(cm20-cm02)*((cm30+cm12)^2-(cm21+cm03)^2)+4*cm11*(cm30+cm12)*(cm21+cm03);  %
    ju(7)=(3*cm21-cm03)*(cm30+cm12)*((cm30+cm12)^2-3*(cm21+cm03)^2)+(cm30-3*cm12)*(cm21+cm03)*(3*(cm30+cm12)^2-(cm21+cm03)^2); 
    qijieju= ju;%abs(log(ju))
    cm = [cm00 cm11 cm02 cm20];
end

三、通过两点计算夹角

python实现

def azimuthAngle( x1,  y1,  x2,  y2):
    angle = 0.0
    dx = x2 - x1
    dy = y2 - y1
    angle = math.atan(dy/dx)
    return (angle * 180 / math.pi)

angle = azimuthAngle(x1,y1,x2,y2)
print("angle=",angle)

 

2015-10-23 19:15:21 lishanlu136 阅读数 6824

学习图像处理,必须首先得理解图像的特征提取,我在这里把图像的特征暂且分为如下五类:面特征、线特征、局部区域特征、点特征和不变点特征。其中,面特征对应图像多分辨率金字塔和图像的矩特征,线特征对应图像的边缘检测,局部区域特征对应图像的斑点特征检测,点特征对应图像的角点检测,不变点特征对应图像的尺度不变特征的提取。下面我一一简单说明一下这些特征。

一、图像金字塔

在数字图像处理领域,多分辨率金字塔化是图像多尺度表示的主要形式。图像金字塔化一般包括两个步骤:图像经过一个低通滤波器进行平滑;然后对这个平滑图像进行抽样,抽样比例一般在水平和垂直方向都是为1/2,从而得到一系列尺寸缩小、分辨率降低的图像。将得到的依次缩小的图像按顺序排列,看上去就像是金字塔,这就是“图像金字塔”的由来。

图像的高斯金字塔:设原图像为G0,以G0作为高斯金字塔的第0层(底层),对原始输入图像进行高斯低通滤波和隔行隔列的降采样,得到高斯金字塔的第1层;再对第1层图像低通滤波和降采样,得到高斯金字塔的第2层;重复以上过程,就可以构成高斯金字塔。可见高斯金字塔的当前层图像是对其前一层图像先进行高斯低通滤波,然后进行隔行隔列的降采样而生成的。当前层图像的大小依次为前一层图像大小的1/4。

图像的拉普拉斯金字塔:它的每一层图像是高斯金字塔本层图像与其高一级的图像(由本层图像下采样得到)经内插放大后图像的差,此过程相当于带通滤波,因此,拉普拉斯金字塔又称为带通金字塔。

MATLAB图像处理工具箱提供了impyramid()函数,用于构造图像的金字塔,其调用格式如下:B=impyramid(A,direction),该函数的功能为:对A进行高斯金字塔变换,direction为'reduce'和'expand',分别对应分解和扩张。

二、图像的矩特征

图像的矩特征是用于可以代表图像的图像描述量,不变矩(Invariant Moments,IMg)是一种高度浓缩的图像特征,具有平移、灰度、尺度、旋转不变性,因此矩和矩函数可被用于图像的模式识别、图像分类、目标识别和场景分析中。

HU矩 几何矩是由Hu(Visual pattern recognition by moment invariants)在1962年提出的,图像f(x,y)的(p+q)阶几何矩定义为 Mpq =∫∫(x^p)*(y^q)f(x,y)dxdy(p,q = 0,1,……∞)矩在统计学中被用来反映随机变量的分布情况,推广到力学中,它被用作刻画空间物体的质量分布。同样的道理,如果我们将图像的灰度值看作是一个二维或三维的密度分布函数,那么矩方法即可用于图像分析领域并用作图像特征的提取。最常用的,物体的零阶矩表示了图像的“质量”:Moo= ∫∫f(x,y )dxdy 一阶矩(M01,M10)用于确定图像质心( Xc,Yc):Xc = M10/M00;Yc = M01/M00;若将坐标原点移至 Xc和 Yc处,就得到了对于图像位移不变的中心矩。如Upq =∫∫[(x-Xc)^p]*[(y-Yc)^q]f(x,y)dxdy。Hu在文中提出了7个几何矩的不变量,这些不变量满足于图像平移、伸缩和旋转不变。如果定义Zpq=Upq/(U20 + U02)^(p+q+2),Hu 的7种矩为:H1=Z20+Z02;H1=(Z20+Z02)^2+4Z11^2;......

三、图像的边缘检测

图像的边缘是指其周围像素灰度急剧变化的那些像素的集合,它是图像最基本的特征。边缘存在于目标、背景和区域中,所以,它是图像分割最重要的依据。边缘检测的基本思想是先检测图像中的边缘点,再按照某种策略将边缘点连接成轮廓,从而构成分割区域。

1、运用一阶微分算子检测图像边缘

一阶微分边缘检测算子也称为梯度边缘算子,它是利用图像在边缘处的阶跃性,即图像梯度在边缘取得极大值的特性进行边缘检测。梯度是一个矢量,它具有方向和模。梯度的模值大小提供了边缘的强度信息,梯度的方向提供了边缘的趋势信息,因为梯度的方向始终是垂直于边缘的方向。在实际使用中,通常利用有限差分进行梯度近似。x方向的导数近似等于I(x+1,y)-I(x,y),y方向的导数近似等于I(x,y+1)-I(x,y),对于3*3模板中心像元的梯度,其梯度可通过下式计算得到:x方向:Mx=(a2+ca3+a4)-(a0+ca7+a6),y方向:My=(a6+ca5+a4)-(a0+ca1+a2),其中a0至a8是由模板左上角a0顺时针旋转至模板中心a8,参数c为加权系数,表示离中心像元较近。当c等于1时,就可以得到Prewitt边缘检测卷积核,当c=2时,就可以得到Sobel边缘检测卷积核。

MATLAB中也提供了相关的图像边缘检测的函数,其调用格式如下:BW = edge(I,'sobel',thresh,direction),BW = edge(I,'prewitt',thresh,direction),

BW = edge(I,'roberts',thresh),其中,I是输入的灰度图像,thresh是阈值,direction是方向。

2、傅里叶变换与图像梯度的关系

实际上,对图像进行二维傅里叶变换得到的频谱图,就是图像梯度的分布图,当然频谱图上的各点与图像上的各点并不存在一一对应的关系,在傅里叶频谱图上看到的明暗不一的亮点,实际上是图像上某一点与领域点差异的强弱,即梯度的大小,也即该点的频率大小(也可以这样理解,图像中的低频部分指低梯度的点,高频部分指高梯度的点),一般来讲,梯度大则该点的亮度强,否则该点亮度弱。通过观察图像的频谱图,我们首先可以看出图像的能量分布,如果频谱图中的暗的点数多,那么实际图像是比较柔和的(因为各点与领域差异不大,梯度相对较小),反之,如果频谱图中亮的点数较多,那么实际图像一定是尖锐的,边界分明且边界两边像素差异较大。

3、运用二阶微分算子检测图像边缘

二阶微分边缘检测算子是利用图像在边缘处的阶跃性导致图像二阶微分在边缘处出现零值这一特性进行边缘检测的,因此,该方法也被称为过零算子和拉普拉斯算子。如下:

x方向二阶微分算子=I(i,j+1)-2I(i,j)+I(i,j-1),y方向二阶微分算子=I(i+1,j)-2I(i,j)+I(i-1,j)

虽然使用二阶微分算子检测边缘的方法简单,但是它是缺点是对噪声十分敏感,同时,也不能提供边缘的方向信息。为了实现对噪声的抑制,Marr等提出了高斯拉普拉斯(LoG)的方法。即采用高斯函数作为低通滤波器对图像滤波后,再对该图像进行二阶微分运算。也可以转换为先对高斯函数进行二阶微分,然后再将结果对图像进行卷积运算。为减少高斯函数二阶微分的计算量,可以直接用高斯差分算子(Difference of Gaussian,DoG)代替它。

在MATLAB中,也提供了相关函数,其调用格式为:BW = edge(I,'log',thresh)。其中,I是输入的灰度图像,thresh是阈值。

4、基于Canny算子检测图像边缘

Canny边缘检测算子是边缘检测算子中最常用的一种,虽然Canny算子也是一阶微分算子,但它对一阶微分算子进行了扩展:主要是在原一阶微分算子的基础上,增加了非最大值抑制和双阈值两项改进。利用非最大值抑制不仅可以抑制多响应边缘,而且还可以提高边缘的定位精度;利用双阈值可以有效地减少边缘的漏检率。Canny算子进行边缘提取主要分为4步:

1)去噪声,通常使用高斯函数对图像进行平滑滤波。

2)计算梯度值与方向角。

3)非最大值抑制

4)滞后阈值化。

在MATLAB中,也提供了相关函数,其调用格式为:BW = edge(I,'canny',thresh)。其中,I是输入的灰度图像,thresh是阈值。

四、斑点特征检测

斑点通常是指与周围有着颜色和灰度差别的区域,斑点检测(Blob Detection)是数字图像处理研究的重要内容,由于斑点代表的是一个区域,相比于单纯的角点,它的稳定性要好,抗噪声能力要强。可作为局部特征中的一个重要特例,斑点在图像配准和立体视觉中充当着重要的角色。

1、利用高斯拉普拉斯(Laplace of Guassian,LoG)算子检测图像斑点是一种十分常用的方法。

注意:图像与某一个二维函数进行卷积运算,实际就是求取图像与这一函数的相似性。

2、DoH斑点

图像像素点的二阶微分Hessian矩阵的行列式DoH(Determinant of Hessian)的值同样也反映了图像局部的结构信息。它对图像中细长结构的斑点具有良好的抑制作用。

第一步:使用不同的δ生成模板,并对图像进行卷积运算;

第二步:在图像的位置空间和尺度空间搜索DoH响应的峰值。

五、角点检测

在图像中,可从两个不同的方向去定义角点:角点是两个边缘的交点;角点是领域内具有两个方向的特征点。角点所在的领域通常也是图像中稳定的、信息丰富的区域。

1、Harris角点

图像Harris角点的检测算法实现步骤归纳为:计算图像I(x,y)在x和y两个方向的梯度Ix、Iy

计算图像两个方向梯度的乘积Ix*Ix,Ix*Iy,Iy*Iy

使用高斯函数对这三个乘积进行高斯加权,生成矩阵M的元素A,B,C,其中M=[A   B

                                                                                                                                                                                                         B   C]                   

计算每个像元的Harris响应值R,并对小于某一阈值的R置零: R = { R: detM- a(tranceM)*(tranceM)< t}

在3*3或5*5的领域内进行非极大值抑制,局部极大值点即为图像中的角点。

在MATLAB中,可以调用C = cornermetric(I,‘Harris’)来检测图像的Harris角点特征,其中,I为输入的灰度图像矩阵;c为角点量度矩阵,用来探测图像中的角点信息,并与I同尺寸,C的值越大表示图像I中的像素越有可能是一个角点。

六、尺度不变特征提取

1、SIFT特征提取

尺度不变特征变换(Scale Invariant Feature Transform,SIFT)是一种图像特征提取与描述算法。SIFT算法由David.G.Lowe于1999年提出并在2004年进行了完善总结。SIFT算法可以处理两幅图像之间发生平移、旋转、尺度变化、光照变化情况下的特征匹配问题,并能在一定程度上对视角变化、仿射变化也具备较为稳定的特征匹配能力。

基本原理与具体实现步骤:

整个过程主要分为四个步骤。

(1)尺度空间峰值选择(Scale space peak selection), 这一步的目的是在尺度空间中选择选择潜在的满足尺度不变性和旋转不变性的关键点(或者称为兴趣点)。

(2)关键点定位(key point localization): 这一步的目的就是精确定位出特征关键点的位置, 涉及到剔除伪关键点。

(3)方向分配(orientation assignment): 这一步的目的就是基于关键点的局部梯度方向,给每个相应的关键点分配方向(Assign orientation to key points)。

(4)关键点描述(Key Point descriptor): 此步骤的目的就是对于每个关键点, 用一个高维度(high dimensional vector,128 维的向量)的向量去描述每个关键点。

                                                                     
2016-07-15 18:17:37 boon_228 阅读数 7304

一 原理

    几何矩是由Hu(Visual pattern recognition by moment invariants)1962年提出的,具有平移、旋转和尺度不变性。 定义如下:

① (p+q)阶不变矩定义

② 对于数字图像,离散化,定义为

 

③ 归一化中心矩定义

④Hu矩定义

 



//#################################################################################//
double M[7] = {0};        //HU不变矩
bool HuMoment(IplImage* img)
{
int bmpWidth = img->width;
int bmpHeight = img->height;
int bmpStep = img->widthStep; 
int bmpChannels = img->nChannels;
    uchar*pBmpBuf = (uchar*)img->imageData;

double m00=0,m11=0,m20=0,m02=0,m30=0,m03=0,m12=0,m21=0;  //中心矩 
double x0=0,y0=0;    //计算中心距时所使用的临时变量(x-x') 
double u20=0,u02=0,u11=0,u30=0,u03=0,u12=0,u21=0;//规范化后的中心矩
//double M[7];    //HU不变矩 
double t1=0,t2=0,t3=0,t4=0,t5=0;//临时变量, 
//double Center_x=0,Center_y=0;//重心 
int Center_x=0,Center_y=0;//重心 
int i,j;            //循环变量

//  获得图像的区域重心(普通矩)
double s10=0,s01=0,s00=0;  //0阶矩和1阶矩  
for(j=0;j<bmpHeight;j++)//y
    {
for(i=0;i<bmpWidth;i++)//x
        {
            s10+=i*pBmpBuf[j*bmpStep+i];
            s01+=j*pBmpBuf[j*bmpStep+i];
            s00+=pBmpBuf[j*bmpStep+i];
        }
    }
    Center_x=(int)(s10/s00+0.5);
    Center_y=(int)(s01/s00+0.5);

//  计算二阶、三阶矩(中心矩)
    m00=s00; 
for(j=0;j<bmpHeight;j++) 
    {
for(i=0;i<bmpWidth;i++)//x 
        { 
            x0=(i-Center_x); 
            y0=(j-Center_y); 
            m11+=x0*y0*pBmpBuf[j*bmpStep+i]; 
            m20+=x0*x0*pBmpBuf[j*bmpStep+i]; 
            m02+=y0*y0*pBmpBuf[j*bmpStep+i]; 
            m03+=y0*y0*y0*pBmpBuf[j*bmpStep+i];
            m30+=x0*x0*x0*pBmpBuf[j*bmpStep+i]; 
            m12+=x0*y0*y0*pBmpBuf[j*bmpStep+i]; 
            m21+=x0*x0*y0*pBmpBuf[j*bmpStep+i]; 
        } 
    } 

//  计算规范化后的中心矩: mij/pow(m00,((i+j+2)/2)
    u20=m20/pow(m00,2); 
    u02=m02/pow(m00,2); 
    u11=m11/pow(m00,2);
    u30=m30/pow(m00,2.5); 
    u03=m03/pow(m00,2.5);
    u12=m12/pow(m00,2.5); 
    u21=m21/pow(m00,2.5);

//  计算中间变量
    t1=(u20-u02); 
    t2=(u30-3*u12); 
    t3=(3*u21-u03); 
    t4=(u30+u12);
    t5=(u21+u03);

//  计算不变矩 
    M[0]=u20+u02; 
    M[1]=t1*t1+4*u11*u11; 
    M[2]=t2*t2+t3*t3; 
    M[3]=t4*t4+t5*t5;
    M[4]=t2*t4*(t4*t4-3*t5*t5)+t3*t5*(3*t4*t4-t5*t5); 
    M[5]=t1*(t4*t4-t5*t5)+4*u11*t4*t5;
    M[6]=t3*t4*(t4*t4-3*t5*t5)-t2*t5*(3*t4*t4-t5*t5);

returntrue;
}

下面的代码计算轮廓的矩,并根据1阶中心矩得到轮廓的质心,代码如下:

src = imread( "../star1.jpg" ,1 );

/// Convert image to gray and blur it
cvtColor( src, src_gray, CV_BGR2GRAY );
blur( src_gray, src_gray, Size(3,3) );

namedWindow( "image", CV_WINDOW_AUTOSIZE );
imshow( "image", src );

Mat canny_output;
vector<vector<Point> > contours;
vector<Vec4i> hierarchy;

//利用canny算法检测边缘
Canny( src_gray, canny_output, thresh, thresh*2, 3 );
namedWindow( "canny", CV_WINDOW_AUTOSIZE );
imshow( "canny", canny_output );
//查找轮廓
findContours( canny_output, contours, hierarchy, CV_RETR_TREE, CV_CHAIN_APPROX_SIMPLE, Point(0, 0) );

//计算轮廓矩
vector<Moments> mu(contours.size() );
for( int i = 0; i < contours.size(); i++ )
    { mu[i] = moments( contours[i], false ); }

//计算轮廓的质心
vector<Point2f> mc( contours.size() );
for( int i = 0; i < contours.size(); i++ )
    { mc[i] = Point2f( mu[i].m10/mu[i].m00 , mu[i].m01/mu[i].m00 ); }

//画轮廓及其质心
Mat drawing = Mat::zeros( canny_output.size(), CV_8UC3 );
for( int i = 0; i< contours.size(); i++ )
    {
    Scalar color = Scalar( rng.uniform(0, 255), rng.uniform(0,255), rng.uniform(0,255) );
    drawContours( drawing, contours, i, color, 2, 8, hierarchy, 0, Point() );
    circle( drawing, mc[i], 4, color, -1, 8, 0 );
    }

namedWindow( "Contours", CV_WINDOW_AUTOSIZE );
imshow( "Contours", drawing );

//打印轮廓面积和轮廓长度
printf("\t Info: Area and Contour Length \n");
for( int i = 0; i< contours.size(); i++ )
    {
    printf(" * Contour[%d] - Area (M_00) = %.2f - Area OpenCV: %.2f - Length: %.2f \n", i, mu[i].m00, contourArea(contours[i]), arcLength( contours[i], true ) );
    Scalar color = Scalar( rng.uniform(0, 255), rng.uniform(0,255), rng.uniform(0,255) );
    drawContours( drawing, contours, i, color, 2, 8, hierarchy, 0, Point() );
    circle( drawing, mc[i], 4, color, -1, 8, 0 );
    }


2014-04-18 13:18:24 cui134 阅读数 2190

利用图像的不变矩特征进行图像匹配是一种基于特征的图像匹配方式,它通过计算两幅图的不变矩特征,然后再计算这两组特征的相似度来判断两幅图是否匹配。

这里利用图像的二阶和三阶矩来导出7个不变矩,这些特征在比例因子小于2和旋转角度不超过45度的条件下,对于平移,旋转,比例因子的变化都是不变的,它反映了图像的固有特性。这7个不变矩的计算公式可以看下面的程序,公式确实写的累啊。

opencv实现的代码如下:

#include "stdafx.h"  
   
#include <opencv2/opencv.hpp>  
#include "highgui.h"  
#include <math.h>  
typedef  unsigned long uint32;
typedef  unsigned int  uint16;
typedef  unsigned char uint8;

IplImage *src_gray1, *src_gray2, *src_gray3;    
IplImage *temp_gray1, *temp_gray2, *temp_gray3;   

CvPoint CalImageCenter( IplImage* I )    //单通道图像,计算图像重心位置
{
	CvPoint CenterPosition;
	double  m10=0;
	double  m00=0;
	double  m01=0;
    int i,j;
	
	for ( i=0; i<I->height; i++ )
	{
		uint8* ptr = (uint8*)( I->imageData + i*I->widthStep );
		for ( j=0; j<I->width; j++ )
		{
			uint8 Pixel = ptr[j];
			m10 += (j)*Pixel;
			m00 += Pixel;
			m01 += (i)*Pixel;		
			
		}
	}
	CenterPosition.x = m10/m00;
	CenterPosition.y = m01/m00;
	return  CenterPosition; 	
}


double CalCenterMoment( IplImage* I, int CenterX, int CenterY, int ip, int jp )  //计算图像的中心矩
{
	int i,j;
	double ICenterMoment=0;
	for ( i=0; i<I->height; i++ )
	{
		uint8* ptr = (uint8*)( I->imageData + i*I->widthStep );
		for ( j=0; j<I->width; j++ )
		{
			uint8 Pixel = ptr[j];
			double temp;       //中间值
			temp = pow( (double)(j-CenterX),ip )*pow( (double)(i-CenterY),jp );
			temp = temp*Pixel;
			ICenterMoment += temp;
		}

	}

	return ICenterMoment;
}

void AllocateImage(IplImage* I,IplImage* T)   //给图像分配大小    
{    
    CvSize sz   = cvGetSize(I);    
    CvSize sz_T = cvGetSize(T);    
        
    
    src_gray1 = cvCreateImage( sz, IPL_DEPTH_8U, 1);    //原图的三个通道    
    src_gray2 = cvCreateImage( sz, IPL_DEPTH_8U, 1);    
    src_gray3 = cvCreateImage( sz, IPL_DEPTH_8U, 1);    
    
    temp_gray1 = cvCreateImage( sz_T, IPL_DEPTH_8U, 1);    //模板的三个通道    
    temp_gray2 = cvCreateImage( sz_T, IPL_DEPTH_8U, 1);    
    temp_gray3 = cvCreateImage( sz_T, IPL_DEPTH_8U, 1);    
    
}    

int _tmain(int argc, _TCHAR* argv[])
{
	IplImage* src_img, *temp_img;                                 //定义变量
	CvPoint   SImageCenter = cvPoint(0,0);
	CvPoint   TImageCenter = cvPoint(0,0);
	double    Su00, Su02, Su20, Su11, Su30, Su12, Su21, Su03;
	double    Tu00, Tu02, Tu20, Tu11, Tu30, Tu12, Tu21, Tu03;
	double    Sa[8], Ta[8];
	double    InterRelateValue=0, S_energy=0, T_energy=0; 
	double    dbr;

	src_img  = cvLoadImage("Images/不变矩.bmp");                 //读入两幅图
	temp_img = cvLoadImage("Images/不变矩模板.bmp");

	AllocateImage( src_img, temp_img);							 //图像初始化

	cvSplit( src_img, src_gray1, src_gray2, src_gray3, 0);        //将两幅三通道图像分解为3幅单通道图像
    cvSplit( temp_img, temp_gray1, temp_gray2, temp_gray3, 0);

	SImageCenter = CalImageCenter(src_gray1);                     //计算两幅图的重心
	TImageCenter = CalImageCenter(temp_gray1);
	
	

	Su00 = CalCenterMoment( src_gray1, 0, 0, 0, 0 );                     //计算两幅图的规格化n阶矩
	Su02 = CalCenterMoment( src_gray1, 0, 0, 0, 2 ) / pow(Su00,2);
	Su20 = CalCenterMoment( src_gray1, 0, 0, 2, 0 ) / pow(Su00,2);
	Su11 = CalCenterMoment( src_gray1, 0, 0, 1, 1 ) / pow(Su00,2);
	Su30 = CalCenterMoment( src_gray1, 0, 0, 3, 0 ) / pow(Su00,2.5);
	Su12 = CalCenterMoment( src_gray1, 0, 0, 1, 2 ) / pow(Su00,2.5);
	Su21 = CalCenterMoment( src_gray1, 0, 0, 2, 1 ) / pow(Su00,2.5);
	Su03 = CalCenterMoment( src_gray1, 0, 0, 0, 3 ) / pow(Su00,2.5);


	Tu00 = CalCenterMoment( temp_gray1, 0, 0, 0, 0 );
	Tu02 = CalCenterMoment( temp_gray1, 0, 0, 0, 2 ) / pow(Tu00,2);
	Tu20 = CalCenterMoment( temp_gray1, 0, 0, 2, 0 ) / pow(Tu00,2);
	Tu11 = CalCenterMoment( temp_gray1, 0, 0, 1, 1 ) / pow(Tu00,2);
	Tu30 = CalCenterMoment( temp_gray1, 0, 0, 3, 0 ) / pow(Tu00,2.5);
	Tu12 = CalCenterMoment( temp_gray1, 0, 0, 1, 2 ) / pow(Tu00,2.5);
	Tu21 = CalCenterMoment( temp_gray1, 0, 0, 2, 1 ) / pow(Tu00,2.5);
	Tu03 = CalCenterMoment( temp_gray1, 0, 0, 0, 3 ) / pow(Tu00,2.5);
	
	Sa[1] = Su20 + Su02;                                                 //计算两幅图的7个不变矩
	Sa[2] = pow( (Su20 - Su02), 2 ) + 4 * pow( Su11, 2 );
	Sa[3] = pow( (Su30 - 3*Su12), 2 ) + pow( (3*Su21-Su03), 2 );
	Sa[4] = pow( (Su30 + Su12), 2 ) + pow( (Su21+Su03), 2 );
	Sa[5] = ( Su30 - 3*Su12 )*( Su30 + Su12 )*( pow( Su30+Su12, 2 ) - 3*pow( Su21+Su03, 2 ) ) 
		    + ( 3*Su21 - Su03 )*( Su21 + Su03 )*( 3*pow( Su30+Su12, 2 ) - pow( Su21+Su03, 2 ) );
	Sa[6] = ( Su20 - Su02 )*( pow( Su30+Su12, 2 ) - pow( Su21+Su03, 2 ) ) + 
			4*Su11*( Su30+Su12 )*( Su21+Su03 );
	Sa[7] = (3*Su21-Su03)*(Su30+Su12)*( pow(Su30+Su12,2) - 3*pow(Su21+Su03,2) ) + 
			(Su30-3*Su12)*(Su21+Su03)*( 3*pow(Su30+Su12,2) - pow(Su21+Su03,2) );

	Ta[1] = Tu20 + Tu02;
	Ta[2] = pow( (Tu20 - Tu02), 2 ) + 4 * pow( Tu11, 2 );
	Ta[3] = pow( (Tu30 - 3*Tu12), 2 ) + pow( (3*Tu21-Tu03), 2 );
	Ta[4] = pow( (Tu30 + Tu12), 2 ) + pow( (Tu21+Tu03), 2 );
	Ta[5] = ( Tu30 - 3*Tu12 )*( Tu30 + Tu12 )*( pow( Tu30+Tu12, 2 ) - 3*pow( Tu21+Tu03, 2 ) ) 
		    + ( 3*Tu21 - Tu03 )*( Tu21 + Tu03 )*( 3*pow( Tu30+Tu12, 2 ) - pow( Tu21+Tu03, 2 ) );
	Ta[6] = ( Tu20 - Tu02 )*( pow( Tu30+Tu12, 2 ) - pow( Tu21+Tu03, 2 ) ) + 
			4*Tu11*( Tu30+Tu12 )*( Tu21+Tu03 );
	Ta[7] = (3*Tu21-Tu03)*(Tu30+Tu12)*( pow(Tu30+Tu12,2) - 3*pow(Tu21+Tu03,2) ) + 
			(Tu30-3*Tu12)*(Tu21+Tu03)*( 3*pow(Tu30+Tu12,2) - pow(Tu21+Tu03,2) );


	for( int i=1;i<8;i++ )
	{
		InterRelateValue += Sa[i]*Ta[i];
		S_energy         += pow(Sa[i],2);
		T_energy         += pow(Ta[i],2);
	}
	dbr = InterRelateValue/sqrt(S_energy*T_energy);                              //计算两幅图7个不变矩的相似度,就是两幅图的相似度
	
	cvNamedWindow("my picture",CV_WINDOW_AUTOSIZE);  
    cvNamedWindow("my model",CV_WINDOW_AUTOSIZE);

	cvShowImage("my picture",src_gray1);             
    cvShowImage("my model",temp_gray1);

	cvWaitKey(0);  
    cvReleaseImage(&src_img);  
    cvReleaseImage(&temp_img);  
    cvReleaseImage(&src_gray1);  
    cvReleaseImage(&src_gray2);  
    cvReleaseImage(&src_gray3);  
    cvReleaseImage(&temp_gray1);  
    cvReleaseImage(&temp_gray2);  
    cvReleaseImage(&temp_gray3); 
	cvDestroyWindow("my picture");  
    cvDestroyWindow("my model");  


	return 0;
}


用于匹配的两幅图如下:

左图的7个不变矩特征Sa[1]~Sa[7]:

右图的7个不变矩特征Ta[1]~Ta[7]:


计算的相似度:



图像矩特征及Hu矩

阅读数 178

图像的矩的知识

阅读数 295

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