本文內容：主要介紹中心矩的定義，以及利用中心矩求取影象主軸方向。

1.4 matlab求取歸一化中心矩、hu不變矩和主軸方向

function test
   
% 	atan(phi) 值域為[-90,90] atan2(y,x) 值域為[-180,180]
% 	2phi範圍為[-180, 180] 因此此處應用atan2
% 	所求主軸方向，為與x軸正向最小夾角，夾角在y正象限，phi>0
% 	y負象限，phi<0
% 	使座標旋轉至水平方向，phi>0時，應沿y負旋轉， phi<0時，應沿
% 	y正旋轉，因此對於旋轉矩陣為：
% 	[cos sin; -sin cos]

    I = imread('BW6.bmp');
  
    [cm ju] = qijieju(uint8(I));
    m00 = cm(1);
    mu11 = cm(2);
    mu02 = cm(3);
    mu20 = cm(4);
    a = mu20 / m00; 
    b = mu11 / m00;
    c = mu02 / m00;
    square = sqrt( 4 * b * b + (a - c) * (a - c) );
	%求主軸方法1
    theta = atan2( 2 * b, a - c + square )*180/pi

    %求主軸方法2
    angle3=atan2(2*mu11,(mu20-mu02))/2*180/pi

end

%求不變矩及中心矩
function [cm ju] = qijieju(I0)
    A=double(I0);
    [nc,nr]=size(A);
    [x,y]=meshgrid(1:nr,1:nc);
    x=x(:);
    y=y(:);
    A=A(:);
    m00=sum(A);
    if m00==0
        m00=eps;
    end
    m10=sum(x.*A);
    m01=sum(y.*A);
    xmean=m10/m00; %重心
    ymean=m01/m00;
    cm00=m00; %歸一化中心矩
    cm02=(sum((y-ymean).^2.*A))/(m00^2);
    cm03=(sum((y-ymean).^3.*A))/(m00^2.5);
    cm11=(sum((x-ymean).*(y-ymean).*A))/(m00^2);
    cm12=(sum((x-ymean).*(y-ymean).^2.*A))/(m00^2.5);
    cm20=(sum((x-xmean).^2.*A))/(m00^2);
    cm21=(sum((x-xmean).^2.*(y-ymean).*A))/(m00^2.5);
    cm30=(sum((x-xmean).^3.*A))/(m00^2.5);
    ju(1)=cm20+cm02;  %
    ju(2)=(cm20-cm02)^2+4*cm11^2; %
    ju(3)=(cm30-3*cm12)^2+(3*cm21-cm03)^2; %
    ju(4)=(cm30+cm12)^2+(cm21+cm03)^2;  %
    ju(5)=(cm30-3*cm12)*(cm30+cm12)*((cm30+cm12)^2-3*(cm21+cm03)^2)+(3*cm21-cm03)*(cm21+cm03)*(3*(cm30+cm12)^2-(cm21+cm03)^2); %
    ju(6)=(cm20-cm02)*((cm30+cm12)^2-(cm21+cm03)^2)+4*cm11*(cm30+cm12)*(cm21+cm03);  %
    ju(7)=(3*cm21-cm03)*(cm30+cm12)*((cm30+cm12)^2-3*(cm21+cm03)^2)+(cm30-3*cm12)*(cm21+cm03)*(3*(cm30+cm12)^2-(cm21+cm03)^2); 
    qijieju= ju;%abs(log(ju))
    cm = [cm00 cm11 cm02 cm20];
end