图像特征中的几何不变矩：Hu矩

对了，我最近开通了微信公众号，计划是两边会同步更新，并逐步的会将博客上的文章同步至公众号中。
感兴趣的朋友可以扫描下方的二维码或者搜索“里先森sements”来关注，欢迎来玩~！

---

本文中的主要参考文献如下：

• 【图像算法】图像特征：几何不变矩–Hu矩 - SkySeraph - 博客园: https://www.cnblogs.com/skyseraph/archive/2011/07/19/2110183.html
• 图像的矩特征 - ☆Ronny丶 - 博客园: https://www.cnblogs.com/ronny/p/3985810.html
• 《图像处理、分析与机器视觉》
• 《数字图像处理》冈萨雷斯–第十一章-不变矩

一、矩的概念

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

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

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

二、Hu矩

2.1 (p+q)阶不变矩的定义

对于灰度分布为$f(x,y)$的图像，其$(p+q)$阶矩定义为：

$$m_{pq}=\iint x^p{y}^{q}f(x,y)dxdyp,q=0,1,2,\cdots$$

其$(p+q)$阶中心矩定义为：

$${\mu }_{pq}=\iint (x-x_0{)}^p(y=y_0{)}^{q}f(x,y)dxdy$$

其中，矩心为：

$$x_0=\frac{m_{10}}{m_{00}},y_0=\frac{m_{01}}{m_{00}}$$

由于数字图像属于离散化数据，离散状态下，$f(x,y)$的$(p+q)$阶普通矩和中心矩的公式如下：

$$m_{pq}=\sum _{x=1}^M\sum _{y=1}^N{x}^p{y}^{q}f(x,y)$$

$${\mu }_{pq}=\sum _{x=1}^M\sum _{y=1}^N(x-x_0{)}^p(y-y_0{)}^{q}f(x,y)\text{其中}p,q=0,1,2,\cdots$$

离散化后的矩心x₀与y₀的计算方法同连续状态下的计算公式。

p、q即在x、y这两个像素位置上的多少次方。在p、q不为0时，求$m_{pq}$可以看作 x的p次方 与 y的q次方 作为权重叠加到当前x、y位置上的灰度值，然后累加。当求$m_{00}$时，可知仅为图像的灰度值进行累加。

2.2 归一化中心距

当图像发生变化，如平移、旋转等，$m_{pq}$也发生变化。而由于${\mu }_{pq}$式与中心有关，其对平移不敏感，但依旧对旋转敏感。为了解决这个问题，这里提出归一化中心矩：

$$y_{pq}=\frac{{\mu }_{pq}}{{\mu }_{00}^r}$$

式中：$$r=\frac{p+q+2}{2}p+q=2,3,\cdots$$

归一化中心矩使用的是中心矩进行的计算。可以注意到，归一化中心矩计算$y_{pq}$时使用的${\mu }_{00}$从数值上来看，就等同于$m_{00}$。

2.3 Hu矩定义

Hu利用二阶和三阶中心矩（即$p+q=2,3$）构造了7个不变矩$I_1\sim I_7$，这7个不变矩构成一组特征量，Hu.M.K在1962年证明了他们具有旋转，缩放和平移不变性,它们在连续图像条件下可以保持平移、缩放和旋转不变。具体定义如下:

$$\{\begin{array}{l}{I}_1=y_{20}+y_{02}\\ \\ I_2=(y_{20}-y_{02}{)}^2+4y_{11}^2\\ \\ I_3=(y_{30}-3y_{12}{)}^2+(3y_{21}-y_{03}{)}^2\\ \\ I_4=(y_{30}+y_{12}{)}^2+(y_{21}+y_{03}{)}^2\\ \\ I_5=(y_{30}-3y_{12})(y_{30}+y_{12})[(y_{30}+y_{12}{)}^2-3(y_{21}+y_{03}{)}^2]+(3y_{21}-y_{03})(y_{21}+y_{03})[3(y_{30}+y_{12}{)}^2-(y_{21}+y_{03}{)}^2]\\ \\ I_6=(y_{20}-y_{02})[(y_{30}+y_{12}{)}^2-(y_{21}+y_{03}{)}^2]+4y_{11}(y_{30}+y_{12})(y_{21}+y_{03})\\ \\ I_7=(3y_{21}-y_{03})(y_{30}+y_{12})[(y_{30}+y_{12}{)}^2-3(y_{21}+y_{03}{)}^2]-(y_{30}-3y_{12})(y_{21}+y_{03})[3(y_{30}+y_{12}{)}^2-(y_{21}+y_{03}{)}^2]\end{array}$$

2.4 Hu矩计算代码

参考自：【图像算法】图像特征：几何不变矩–Hu矩 - SkySeraph - 博客园，修正了部分格式问题与变量名，方便对照上面进行查看。

double I[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 u00=0,u11=0,u20=0,u02=0,u30=0,u03=0,u12=0,u21=0;     //中心矩 
    double x0=0,y0=0;                                           //计算中心距时所使用的临时变量（x-x'） 
    double y20=0,y02=0,y11=0,y30=0,y03=0,y12=0,y21=0;           //规范化后的中心矩
    double t1=0,t2=0,t3=0,t4=0,t5=0;                            //临时变量， 
    int center_x0=0,center_y0=0;                                //矩心 
    int i,j;                                                    //循环变量

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

    //  计算二阶、三阶矩(中心矩)
    u00=m00;                        //归一化中心矩计算时使用的u00从数值上来看，就等同于m00。
    for(j=0;j<bmpHeight;j++) 
    {
        for(i=0;i<bmpWidth;i++)//x 
        { 
            x0=(i-center_x0); 
            y0=(j-center_x0); 
            u11+=x0*y0*pBmpBuf[j*bmpStep+i]; 
            u20+=x0*x0*pBmpBuf[j*bmpStep+i]; 
            u02+=y0*y0*pBmpBuf[j*bmpStep+i]; 
            u03+=y0*y0*y0*pBmpBuf[j*bmpStep+i];
            u30+=x0*x0*x0*pBmpBuf[j*bmpStep+i]; 
            u12+=x0*y0*y0*pBmpBuf[j*bmpStep+i]; 
            u21+=x0*x0*y0*pBmpBuf[j*bmpStep+i]; 
        } 
    } 

    //  计算规范化后的中心矩: uij/pow(u00,((i+j+2)/2) , 而u00 = m00，故直接用m00
    y20=u20/pow(m00,2); 
    y02=u02/pow(m00,2); 
    y11=u11/pow(m00,2);
    y30=u30/pow(m00,2.5); 
    y03=u03/pow(m00,2.5);
    y12=u12/pow(m00,2.5); 
    y21=u21/pow(m00,2.5);

    //  计算中间变量
    t1=(y20-y02); 
    t2=(y30-3*y12); 
    t3=(3*y21-y03); 
    t4=(y30+y12);
    t5=(y21+y03);

    //  计算不变矩 
    I[0]=y20+y02; 
    I[1]=t1*t1+4*y11*y11; 
    I[2]=t2*t2+t3*t3; 
    I[3]=t4*t4+t5*t5;
    I[4]=t2*t4*(t4*t4-3*t5*t5)+t3*t5*(3*t4*t4-t5*t5); 
    I[5]=t1*(t4*t4-t5*t5)+4*y11*t4*t5;
    I[6]=t3*t4*(t4*t4-3*t5*t5)-t2*t5*(3*t4*t4-t5*t5);

    return true;
}