图像矩

https://www.cnblogs.com/ronny/p/3985810.html

1. 矩的概念

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

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

比较重要的有两种情况：

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

2. c=E(X)c=E(X)。这时μk=E[(X−EX)k]μk=E[(X−EX)k]称为XX的kk阶中心矩。

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

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

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

2. Hu矩

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

 

mpq=∑i=1M∑j=1Nipjqf(i,j)mpq=∑i=1M∑j=1Nipjqf(i,j)

 

 

μpq=∑i=1M∑j=1N(i−i¯)p(j−j¯)qf(i,j)μpq=∑i=1M∑j=1N(i−i¯)p(j−j¯)qf(i,j)

 

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

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

μ00=∑Mi=1∑Nj=1(i−i¯)0(j−j¯)0f(i,j)=m00μ00=∑i=1M∑j=1N(i−i¯)0(j−j¯)0f(i,j)=m00

μ10=∑Mi=1∑Nj=1(i−i¯)1(j−j¯)0f(i,j)=0μ10=∑i=1M∑j=1N(i−i¯)1(j−j¯)0f(i,j)=0

μ01=∑Mi=1∑Nj=1(i−i¯)0(j−j¯)1f(i,j)=0μ01=∑i=1M∑j=1N(i−i¯)0(j−j¯)1f(i,j)=0

μ11=∑Mi=1∑Nj=1(i−i¯)1(j−j¯)1f(i,j)=m11−y¯m10μ11=∑i=1M∑j=1N(i−i¯)1(j−j¯)1f(i,j)=m11−y¯m10

μ20=∑Mi=1∑Nj=1(i−i¯)2(j−j¯)0f(i,j)=m20−y¯m01μ20=∑i=1M∑j=1N(i−i¯)2(j−j¯)0f(i,j)=m20−y¯m01

μ02=∑Mi=1∑Nj=1(i−i¯)0(j−j¯)2f(i,j)=m02−y¯m01μ02=∑i=1M∑j=1N(i−i¯)0(j−j¯)2f(i,j)=m02−y¯m01

μ30=∑Mi=1∑Nj=1(i−i¯)3(j−j¯)0f(i,j)=m30−2x¯m20+2x¯2m10μ30=∑i=1M∑j=1N(i−i¯)3(j−j¯)0f(i,j)=m30−2x¯m20+2x¯2m10

μ12=∑Mi=1∑Nj=1(i−i¯)1(j−j¯)2f(i,j)=m12−2y¯m11−x¯m02+2y¯2m10μ12=∑i=1M∑j=1N(i−i¯)1(j−j¯)2f(i,j)=m12−2y¯m11−x¯m02+2y¯2m10

μ21=∑Mi=1∑Nj=1(i−i¯)2(j−j¯)1f(i,j)=m21−2x¯m11−y¯m20+2x¯2m01μ21=∑i=1M∑j=1N(i−i¯)2(j−j¯)1f(i,j)=m21−2x¯m11−y¯m20+2x¯2m01

μ03=∑Mi=1∑Nj=1(i−i¯)0(j−j¯)3f(i,j)=m03−2y¯m02+2y¯2m01μ03=∑i=1M∑j=1N(i−i¯)0(j−j¯)3f(i,j)=m03−2y¯m02+2y¯2m01

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

 

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

 

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

Φ1=η20+η02Φ1=η20+η02

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

Φ3=(η20−3η12)2+3(η21−η03)2Φ3=(η20−3η12)2+3(η21−η03)2

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

Φ5=(η30+3η12)(η30+η12)[(η30+η12)2−3(η21+η03)2]+(3η21−η03)(η21+η03)[3(η30+η12)2−(η21+η03)2]Φ5=(η30+3η12)(η30+η12)[(η30+η12)2−3(η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)Φ6=(η20−η02)[(η30+η12)2−(η21+η03)2]+4η11(η30+η12)(η21+η03)

Φ7=(3η21−η03)(η30+η12)[(η30+η12)2−3(η21+η03)2]+]+(3η12−η30)(η21+η03)[3(η30+η12)2−(η21+η03)2]Φ7=(3η21−η03)(η30+η12)[(η30+η12)2−3(η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|log|Φi|

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

   

类别	log|Φ1|log|Φ1|	log|Φ2|log|Φ2|	log|Φ3|log|Φ3|	log|Φ4|log|Φ4|	log|Φ5|log|Φ5|	log|Φ6|log|Φ6|	log|Φ7|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)W(x)在区间(a,b)(a,b)可积，且W(x)≥0W(x)≥0，则可作为权函数。

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

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

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

 

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

 

Rnm(ρ)=∑s=0(n−|m|)/2(−1)s[(n−s)!]ρn−2ss!(n+|m|2−s)!(n+|m|2+s)!Rnm(ρ)=∑s=0(n−|m|)/2(−1)s[(n−s)!]ρn−2ss!(n+|m|2−s)!(n+|m|2+s)!

 

 

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

 

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

Zernike多项式Vnm(x,y)=Vnm(ρ,θ)Vnm(x,y)=Vnm(ρ,θ)是定义在单位圆x2+y2≤1x2+y2≤1上的正交复函数的集合，具有重要的递推性质，即RnmRnm可由R(n−2)mR(n−2)m和R(n−4)mR(n−4)m得到，公式如下：

 

Rnm(ρ)=[(K22ρ2+K3)R(n−2)m(ρ)+K4R(n−4)m(ρ)]K1Rnm(ρ)=[(K22ρ2+K3)R(n−2)m(ρ)+K4R(n−4)m(ρ)]K1

 

 

Rmm(ρ)=ρmRmm(ρ)=ρm

 

式中：K1=(n+1)(n−1)(n−2)/2,K2=2n(n−1)(n−2),K3=−(n−1)3,K4=−n(n−1)(n−3)/2K1=(n+1)(n−1)(n−2)/2,K2=2n(n−1)(n−2),K3=−(n−1)3,K4=−n(n−1)(n−3)/2。

4.2 Zernike矩的定义

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

 

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

 

上式中的ZnmZnm就是Zernike矩。

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

 

Znm=n+1π∫10∫2π0[Vnm(ρ,θ)]f(ρ,θ)ρdydxdρdθZnm=n+1π∫01∫02π[Vnm(ρ,θ)]f(ρ,θ)ρdydxdρdθ

 

 

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

 

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

4.3 Zernike矩的计算

从Zernike矩的计算公式上来看，对于二维图像，其Zernike矩ZnmZnm为复数，将其实部和虚部分别记为CnmCnm和SnmSnm，则有：

 

Cnm=2n+2π∫10∫2π0[Rnm(ρ)cos(mθ)f(ρ,θ)ρdρdθCnm=2n+2π∫01∫02π[Rnm(ρ)cos(mθ)f(ρ,θ)ρdρdθ

 

 

Cnm=2n+2π∫10∫2π0[Rnm(ρ)sin(mθ)f(ρ,θ)ρdρdθCnm=2n+2π∫01∫02π[Rnm(ρ)sin(mθ)f(ρ,θ)ρdρdθ

 

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

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

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

如果|x|=r|x|=r，则：

 

σ=2(r−x)y|y|+xyrσ=2(r−x)y|y|+xyr

 

如果|y|=r|y|=r，则：

 

σ=2y−xyrσ=2y−xyr

 

我们容易计算出，rr的取值范围为1∼N/21∼N/2，σσ的取值范围是1∼8r1∼8r，再根据参数(r,σ)(r,σ)定义相应的极坐标：

 

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

 

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

 

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

 

 

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

 

1. 确定图像的大小N×NN×N，即公式中的NN；

2. 确定r,σr,σ的范围；

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

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

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

类别	log|Z11|log|Z11|	log|Z20|log|Z20|	log|Z22|log|Z22|	log|Z31|log|Z31|	log|Z40|log|Z40|	log|Z42|log|Z42|	log|Z44|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. 计算待识别对象的不变矩，输入神经网络就可得到待识别对象的类型，或者计算待识别对象不变矩与类别对象不变矩之间的距离，选择最短距离的类别作为待识别对象的类别。

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

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

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