moment1.m

使用matlab软件编程

clear
% % % % 矩的计算
code=imread('0013.bmp');
% figure,imshow(code);
[high1 wide1]=size(code);
% % % % % % % % % % % % % % % % % % 
% % % % % % % % % % % % % % % % % % 
ndots=0;
xav=0;
yav=0;
for x=1:high1
    for y=1:wide1
        if code(x,y)==1
            ndots=ndots+1;
            xav=xav+x;
            yav=yav+y;
        end

    end
end
xav=xav/ndots;
yav=yav/ndots;
% % % % % % % % % % % % % % % % % % 求出x，y的均值
% % % % % % % % % % % % % % % % % % 
% % % add u02,u20
u02=0;u20=0;u11=0;u30a=0;u30s=0;u03a=0;u03s=0;
u12a=0;u12s=0;u21a=0;u21s=0;
xmax=0;xmin=0;ymax=0;ymin=0;

for x=1:high1
    for y=1:wide1
        if code(x,y)==1
            if x>xmax
                xmax=x;
            end
            if x<xmin
                xmin=x;
            end
            if y>ymax
                ymax=y;
            end
            if y<ymin
                ymin=y;
            end
u02=u02+(y-yav)^2;
u20=u20+(x-xav)^2;
u11=u11+(y-yav)*(x-xav);
u30=(x-xav)^3;
if u30>0
    u30a=u30a+u30;
else
    u30s=u30s-u30;
end

u03=(y-yav)^3;
    if u03>0
        u03a=u03a+u03;
    else
        u03s=u03s-u03;
    end

 u12=(x-xav)*(y-yav)*(y-yav);
    if u12>0
        u12a=u12a+u12;
    else
        u12s=u12s-u12;
    end
    
    u21=(x-xav)*(x-xav)*(y-yav);
    if u21>0
        u21a=u21a+u21;
    else
        u21s=u21s-u21;
    end
        end
    end
end
   
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
%%%以上计算出了各阶矩
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
theta=atan((u02-u20+((u20-u02)*(u20-u02)+4*u11*u11)^0.5)/2*u11);%%%%计算θ；
lm1=((u20+u02)+((u20-u02)*(u20-u02)+4*u11*u11)^0.5)/2;
lm2=((u20+u02)-((u20-u02)*(u20-u02)+4*u11*u11)^0.5)/2;%%%%%%%%%%%%%%计算λ1λ2；
size=((xmax-xmin)*(ymax-ymin))^0.5;


fg1=theta*2/pi;%%%%orientation
fg2=(lm1-lm2)/(lm1+lm2);%%%%inertial ratio
fg3=((u20-u02)/(u20+u02)+1)/2;%%%%aspect ration
fg4=((u20+u02)/ndots)^0.5/size;%%%%%%spreadness

fg5=((u30a-u30s)/(u30a+u30s)+1)/2;
fg6=((u03a-u03s)/(u03a+u03s)+1)/2;
fg7=((u21a-u21s)/(u21a+u21s)+1)/2;
fg8=((u12a-u12s)/(u12a+u12s)+1)/2;
% fg=[fg1 fg2 fg3 fg3 fg5 fg6 fg7 fg8]
u=[u11 u20 u02 u12 u21 u03 u30]