Hu矩

close all;
clear all;
I1=imread('lena.bmp');
angle=30;
T=[cos(angle),sin(angle),0;-sin(angle),cos(angle),0;0,0,1];
t1 = maketform('affine',T);
I2 = imtransform(I1,t1,'bilinear','fill',0);%旋转
I2=imresize(I2,0.5);%缩放
figure;
imshow(I1);
figure;
imshow(I2);
H1=Humoment(I1)
H2=Humoment(I2)

function Hu = Humoment(I)
I=double(I);
m00=sum(sum(I));
m10=0;m01=0;
[row,col]=size(I);
for i=1:row
    for j=1:col
        m10=m10+i*I(i,j);
        m01=m01+j*I(i,j);
    end
end
xb=m10/m00;%重心坐标
yb=m01/m00;

%中心距
u02 = 0;u03 = 0;u11 = 0;u12 = 0;u20 = 0;u21 = 0;u30 = 0;
for i=1:row
    for j=1:col        
        u02=u02+(j-yb)^2*I(i,j);
        u03=u03+(j-yb)^3*I(i,j);
        u11=u11+(i-xb)*(j-yb)*I(i,j);
        u12=u12+(i-xb)*(j-yb)^2*I(i,j);
        u20=u20+(i-xb)^2*I(i,j);
        u21=u21+(i-xb)^2*(j-yb)*I(i,j);
        u30=u30+(i-xb)^3*I(i,j);
    end
end

%归一化中心距
n02=u02/m00^2;
n03=u03/m00^2.5;
n11=u11/m00^2;
n12=u12/m00^2.5;
n20=u20/m00^2;
n21=u21/m00^2.5;
n30=u30/m00^2.5;        

%Hu矩
h1 = n20 + n02;                     
h2 = (n20-n02)^2 + 4*(n11)^2;
h3 = (n30-3*n12)^2 + (3*n21-n03)^2;  
h4 = (n30+n12)^2 + (n21+n03)^2;
h5 = (n30-3*n12)*(n30+n12)*((n30+n12)^2-3*(n21+n03)^2)+(3*n21-n03)*(n21+n03)*(3*(n30+n12)^2-(n21+n03)^2);
h6 = (n20-n02)*((n30+n12)^2-(n21+n03)^2)+4*n11*(n30+n12)*(n21+n03);
h7 = (3*n21-n03)*(n30+n12)*((n30+n12)^2-3*(n21+n03)^2)-(3*n12-n30)*(n21+n03)*(3*(n30+n12)^2-(n21+n03)^2);
Hu = [h1 h2 h3 h4 h5 h6 h7];   

结果：

对于两图像做测试，其中第二幅是第一幅旋转30度并缩小到1/4大小：

上述两幅图计算出的Hu矩都为：

 0.0012    0.0000    0.0000    0.0000   -0.0000   -0.0000    0.0000