图像的矩特征

1. 矩的概念

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

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

比较重要的有两种情况：

1. $c=0$。这时$a_k=E(X^k)$称为$X$的$k$阶原点矩

2. $c=E(X)$。这时$\mu_k=E[(X-EX)^k]$称为$X$的$k$阶中心矩。

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

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

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

2. Hu矩

一幅$M\times N$的数字图像$f(i,j)$，其$p+q$阶几何矩$m_{pq}$和中心矩$\mu_{pq}$为：

$$m_{pq}=\sum_{i=1}^M\sum_{j=1}^Ni^pj^qf(i,j)$$

$$\mu_{pq}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^p(j-\bar{j})^qf(i,j)$$

其中$f(i,j)$为图像在坐标点$(i,j)$处的灰度值。$\bar{i}=m_{10}/m_{00},\bar{j}=m_{01}/m_{00}$

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

$\mu{00}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^0(j-\bar{j})^0f(i,j)=m_{00}$

$\mu{10}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^1(j-\bar{j})^0f(i,j)=0$

$\mu{01}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^0(j-\bar{j})^1f(i,j)=0$

$\mu{11}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^1(j-\bar{j})^1f(i,j)=m_{11}-\bar{y}m_{10}$

$\mu{20}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^2(j-\bar{j})^0f(i,j)=m_{20}-\bar{y}m_{01}$

$\mu{02}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^0(j-\bar{j})^2f(i,j)=m_{02}-\bar{y}m_{01}$

$\mu{30}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^3(j-\bar{j})^0f(i,j)=m_{30}-2\bar{x}m_{20}+2\bar{x}^2m_{10}$

$\mu{12}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^1(j-\bar{j})^2f(i,j)=m_{12}-2\bar{y}m_{11}-\bar{x}m_{02}+2\bar{y}^2m_{10}$

$\mu{21}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^2(j-\bar{j})^1f(i,j)=m_{21}-2\bar{x}m_{11}-\bar{y}m_{20}+2\bar{x}^2m_{01}$

$\mu{03}=\sum_{i=1}^M\sum_{j=1}^N(i-\bar{i})^0(j-\bar{j})^3f(i,j)=m_{03}-2\bar{y}m_{02}+2\bar{y}^2m_{01}$

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

$$\eta_{pq}=\frac{\mu_{pa}}{\mu_{00}^\gamma},(\gamma=\frac{p+q}{2},p+q=2,3,\dots)$$

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

$\Phi_1=\eta_{20}+\eta_{02}$

$\Phi_2=(\eta_{20}-\eta_{02})^2+4\eta_{11}^2$

$\Phi_3=(\eta_{20}-3\eta_{12})^2+3(\eta_{21}-\eta_{03})^2$

$\Phi_4=(\eta_{30}+\eta_{12})^2+(\eta_{21}+\eta_{03})^2$

$\Phi_5=(\eta_{30}+3\eta_{12})(\eta_{30}+\eta_{12})[(\eta_{30}+\eta_{12})^2-3(\eta_{21}+\eta_{03})^2]+(3\eta_{21}-\eta_{03})(\eta_{21}+\eta_{03})[3(\eta_{30}+\eta_{12})^2-(\eta_{21}+\eta_{03})^2]$

$\Phi_6=(\eta_{20}-\eta_{02})[(\eta_{30}+\eta_{12})^2-(\eta_{21}+\eta_{03})^2]+4\eta_{11}(\eta_{30}+\eta_{12})(\eta_{21}+\eta_{03})$

$\Phi_7=(3\eta_{21}-\eta_{03})(\eta_{30}+\eta_{12})[(\eta_{30}+\eta_{12})^2-3(\eta_{21}+\eta_{03})^2]+]+(3\eta_{12}-\eta_{30})(\eta_{21}+\eta_{03})[3(\eta_{30}+\eta_{12})^2-(\eta_{21}+\eta_{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|\Phi_i|$

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

类别	$log|\Phi_1|$	$log|\Phi_2|$	$log|\Phi_3|$	$log|\Phi_4|$	$log|\Phi_5|$	$log|\Phi_6|$	$log|\Phi_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)\ge 0$，则可作为权函数。

对于一个多项式的序列$f_i$和权函数$W(x)$，定义内积：$<f_m,f_n>=\int_{a}^{b}f_m(x)f_n(x)W(x)dx$

若$n \ne m,<f_m,f_n>=0$，这些多项式则称为正交多项式。若$f_i$除了正交之外，更有$<f_m,f_n>=1$的话，则称为规范正交多项式。

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

 

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

$$R_{nm}(\rho) = \sum_{s=0}^{(n-|m|)/2}\frac{(-1)^s[(n-s)!]\rho^{n-2s}}{s!(\frac{n+|m|}{2}-s)!(\frac{n+|m|}{2}+s)!}$$

$$V_{nm}(x,y)=V_{nm}(\rho,\theta)=R_{nm}(\rho)e^{jm\theta}$$

其中$R_{nm}(\rho)$表示点$(x,y)$的径向多项式，$V_{nm}(x,y)$为Zernike正交多项式，$n,m$为正交多项式的阶数，$n$是非负整数，$n-|m|$是偶数，并且$n\ge|m|$。

Zernike多项式$V_{nm}(x,y)=V_{nm}(\rho,\theta)$是定义在单位圆$x^2+y^2\le 1$上的正交复函数的集合，具有重要的递推性质，即$R_{nm}$可由$R_{(n-2)m}$和$R_{(n-4)m}$得到，公式如下：

$$R_{nm}(\rho)=\frac{[(K_2^2\rho^2+K_3)R_{(n-2)m}(\rho)+K_4R_{(n-4)m}(\rho)]}{K_1}$$

$$R_{mm}(\rho)=\rho^m$$

式中：$K_1=(n+1)(n-1)(n-2)/2,K_2=2n(n-1)(n-2),K_3=-(n-1)^3,K_4=-n(n-1)(n-3)/2$。

4.2 Zernike矩的定义

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

$$f(x,y)=\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}Z_{nm}V_{n,m}(\rho,\theta)$$

上式中的$Z_{nm}$就是Zernike矩。

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

$$Z_{nm}=\frac{n+1}{\pi}\int_{0}^{1}\int_{0}^{2\pi}[V_{nm}(\rho,\theta)]f(\rho,\theta)\rho\frac{dy}{dx}d\rho d\theta$$

$$=\frac{n+1}{\pi}\iint R_{nm}(\rho)e^{jm\theta}f(\rho,\theta)d\rho d\theta$$

式中$\rho=\sqrt{x^2+y^2}(-1<x,y<1)$，$\theta$为轴$x$与$\rho$矢量在逆时针方向的夹角；$R_{nm}(\rho)$表示点$(x,y)$的径向多项式。

4.3 Zernike矩的计算

从Zernike矩的计算公式上来看，对于二维图像，其Zernike矩$Z_{nm}$为复数，将其实部和虚部分别记为$C_{nm}$和$S_{nm}$，则有：

$$C_{nm}=\frac{2n+2}{\pi}\int_{0}^{1}\int_{0}^{2\pi}[R_{nm}(\rho)cos(m\theta)f(\rho,\theta)\rho d\rho d\theta$$

$$C_{nm}=\frac{2n+2}{\pi}\int_{0}^{1}\int_{0}^{2\pi}[R_{nm}(\rho)sin(m\theta)f(\rho,\theta)\rho d\rho d\theta$$

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

对于$N\times N$的图像$f(x,y)$，令坐标原点位于图像的中心，则$-N/2\le x,y\le N/2$，对于像素$(x,y)$，引入2个参数$(r,\sigma)$，唯一对应于像素，其定义为：

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

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

$$\sigma=\frac{2(r-x)y}{|y|}+\frac{xy}{r}$$

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

$$\sigma=2y-\frac{xy}{r}$$

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

$$\rho=2r/N,\theta=\pi\sigma(4r)$$

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

$$C_{nm}=\frac{2n+2}{N^2}\sum_{r=1}^{N/2}R_{nm}(2r/N)\sum_{\sigma=1}^{8r}cos\frac{\pi m\sigma}{4r}f(r,\sigma)$$

$$S_{nm}=\frac{2n+2}{N^2}\sum_{r=1}^{N/2}R_{nm}(2r/N)\sum_{\sigma=1}^{8r}sin\frac{\pi m\sigma}{4r}f(r,\sigma)$$

1. 确定图像的大小$N\times N$，即公式中的$N$；

2. 确定$r,\sigma$的范围；

3. 利用Zernike多项式的递推性质计算各阶$R_{nm}(\rho)$，并结合上面Zernike矩计算公式，算出$C_{nm},S_{nm}$

4. 对$C_{nm},S_{nm}$求模，进而计算得到$|Z_{nm}|$

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

类别	$log|Z_{11}|$	$log|Z_{20}|$	$log|Z_{22}|$	$log|Z_{31}|$	$log|Z_{40}|$	$log|Z_{42}|$	$log|Z_{44}|$
原图	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++版的实现补上。

function        [A_nm,zmlist,cidx,V_nm]        = zernike(img,n,m) 
% 功能：计算输入图像的Zernike矩 
% 输入：img-灰度图像 
%       n-阶数 
% 输出：V_nm-n阶的Zernike多项式，定义为在极坐标系中p，theta的函数 
%       cidx-表示虚部值 
%       A_nm-Zernike矩 

if nargin>0 
if nargin==1 
        n=0; 
end 
d=size(img);                        
img=double(img); 
% 取步长 
xstep=2/(d(1)-1); 
ystep=2/(d(2)-1); 
% 画方格 
[x,y]=meshgrid(-1:xstep:1,-1:ystep:1); 
circle1= x.^2 + y.^2; 
% 提取符合circle1<=1的数 
inside=find(circle1<=1); 
% 构造size（d）*size(d)的矩阵 
mask=zeros(d); 
% 构造size(inside)*size(inside)的全为1的矩阵赋值给mask（inside） 
mask(inside)=ones(size(inside)); 
% 计算图像的复数表示 
[cimg,cidx]=clipimg(img,mask); 
% 计算Z的实部和虚部 
z=clipimg(x+i*y,mask); 
% 计算复数的模，sqrt(x,y),z=x+iy 
p=0.9*abs(z);   ; 
% 计算复数z的辐角值（tanz） 
theta=angle(z); 
c=1; 
for order=1:length(n) 
        n1=n(order); 
           if nargin<3 
                m=zpossible(n1);                
        end 
        for r=1:length(m) 
  V_nmt=zpoly(n1,m(r),p,theta); 
% conj是求复数的共轭 
zprod=cimg.*conj(V_nmt); 
% (n1+1)/π*∑∑(zprod); 对于图像而言求和代替了求积分 
A_nm(c)=(n1+1)*sum(sum(zprod))/pi; 
                zmlist(c,1:2)=[n1 m(r)]; 
                if nargout==4 
                        V_nm(:,c)=V_nmt; 
                end 
                c=c+1; 
        end 
end 
else 
end 

%%%%%子函数%%%%% 
function [cimg,cindex,dim]=clipimg(img,mask) 
%功能：计算复数的实部和虚部 
dim=size(img);  
cindex=find(mask~=0); 
cimg=img(cindex); 
return; 

%%%%%子函数%%%%% 
function  [m]=zpossible(n) 
% 功能：判断n是偶数还是奇数，是偶数时，m取0,2,4,6等,否则取奇数赋值m 
if iseven(n) 
        m=0:2:n; 
else 
        m=1:2:n; 
end 
return; 

%%%%%子函数%%%%% 
function        [V_nm,mag,phase]=zpoly(n,m,p,theta) 
%功能：计算Zernike矩多项式 
R_nm=zeros(size(p)); % 产生size(p)*size(p)的零矩阵赋给R_nm 
a=(n+abs(m))/2; 
b=(n-abs(m))/2; 
total=b; 
for s=0:total 
        num=((-1)^s)*fac(n-s)*(p.^(n-2*s)); % (-1).-1*(n-s)!r.^(n-2*s) 
        den=fac(s)*fac(a-s)*fac(b-s); % s!*(a-s)!*(b-s)! 
        R_nm=R_nm + num/den; % R_nm是一个实数值的径向多项式 
end 
mag=R_nm; % 赋值 
phase=m*theta; 
V_nm=mag.*exp(i*phase); % V_nm为n阶的Zernike多项式，定义为在极坐标系中p，theta的函数 
return; 

%%%%%子函数%%%%% 
function [factorial]=fac(n) 
%功能：求n的阶乘 
maxno=max(max(n)); 
zerosi=find(n<=0); %取n小于等于0的数 
n(zerosi)=ones(size(zerosi)); 
factorial=n; 
findex=n; 
for i=maxno:-1:2 
        cand=find(findex>2); 
        candidates=findex(cand); 
        findex(cand)=candidates-1; 
        factorial(cand)=factorial(cand).*findex(cand); 
end 
return; 

function [verdict]=iseven(candy) 
verdict=zeros(size(candy)); 
isint=find(isint(candy)==1); 
divided2=candy(isint)/2; 
evens=(divided2==floor(divided2)); 
verdict(isint)=evens; 
return; 

function [verdict]=isint(candy) 
verdict        = double(round(candy))==candy; 
return;

5. 总结

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

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

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

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

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

6. 参考资料

[1] 《现代数字图像处理》（matlab版）

[2] 正交多项式WIKI

[3] opencv形态描述