图像处理基础知识（二）—— 中心矩求主轴方向

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