图像的矩特征

2014-09-23 15:33  ☆Ronny丶  阅读(42266)  评论(18编辑  收藏  举报

1. 矩的概念

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

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

2. Hu矩

$$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)$$

$\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)$$

$\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;
}

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

参数说明：

• 输入参数：moments即为上面一个函数计算得到的moments类型。
• 输出参数：hu是一个含有7个数的数组。
int main(int argc, char** argv)
{
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_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矩多项式

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}$$

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$$

4.2 Zernike矩的定义

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

$$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$$

4.3 Zernike矩的计算

$$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$$

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

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

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

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

$$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}|$

 类别 $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

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)的矩阵
% 计算图像的复数表示
% 计算Z的实部和虚部
% 计算复数的模，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

%%%%%子函数%%%%%
%功能：计算复数的实部和虚部
dim=size(img);
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;
View Code

5. 总结

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

6. 参考资料

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

[2] 正交多项式WIKI

[3] opencv形态描述