kblmi.m

图像对齐（inage registration）工具箱

function b=kblmi(a,r)
% Computing complex moment invariants (up to r-th order, r < 8)
% for f*h, h(x,y) symmetric with respect to both axes (AS blur)
% 
% b = kblmi(a,r)
%
% a(n,n) - image matrix
% r - r-th order

vel=size(a);
mc=kmr(a,r);   %spocte jednotlive komplexni momenty
% mc is a matrix of complex moments up to the p-th order of the image img
% \c{pq} = mc(p+1,q+1)

mc00 =mc(1,1);
mc002=mc00*mc00;
mc003=mc002*mc00;

mc10=mc(2,1);
mc01=mc(1,2);

b1=[mc10 mc01];

pq=1;

b1=[mc10 mc01]/((vel(1)/2)^(pq));



% 1 2 
b2=[];
b3=[];
b4=[];
b5=[];
b6=[];

if r>2
	mc11=mc(2,2);
	mc21=mc(3,2);
	mc12=mc(2,3);
	mc30=mc(4,1);
	mc03=mc(1,4);
	mc20=mc(3,1);
	mc02=mc(1,3);

	mc312=mc12 -   (2*mc01*mc11+mc10*mc02)/mc00;

	mc321=mc21 -   (2*mc10*mc11+mc01*mc20)/mc00;

	mc330=mc30 - 3*(mc10*mc20)/mc00;

	mc303=mc03 - 3*(mc01*mc02)/mc00;


b2=[mc303 mc330 mc312 mc321 ];

pq=3;
b2=[mc303 mc330 mc312 mc321 ]/((vel(1)/2)^(pq));
% 3 4 5 6
if r>4
	mc50=mc(6,1);
	mc05=mc(1,6);
	mc40=mc(5,1);
	mc04=mc(1,5);
	mc41=mc(5,2);
	mc14=mc(2,5);
	mc32=mc(4,3);
	mc23=mc(3,4);
	mc31=mc(4,2);
	mc13=mc(2,4);
	mc22=mc(3,3);

	mc202=mc20*mc20;
	mc022=mc02*mc02;
	mc112=mc11*mc11;

	mc550=mc50 - 5*(2*mc20*mc30+mc10*mc40)/mc00+30*(mc10*mc202)/mc002;

	mc505=mc05 - 5*(2*mc02*mc03+mc01*mc04)/mc00+30*(mc01*mc022)/mc002;


	mc541=mc41 - (6*mc20*mc21+mc01*mc40+4*mc11*mc30+4*mc31*mc10)/mc00 + 6*(4*mc10*mc11*mc20+mc202*mc01)/mc002;

	mc514=mc14 - (6*mc02*mc12+mc10*mc04+4*mc11*mc03+4*mc13*mc01)/mc00 + 6*(4*mc01*mc11*mc02+mc022*mc10)/mc002;


	mc532=mc32 - 	(3*mc20*mc12+6*mc11*mc21+2*mc01*mc31+mc02*mc30+3*mc10*mc22)/mc00 + 6*(2*mc20*mc01*mc11+mc20*mc10*mc02+2*mc10*mc112)/mc002;

	mc523=mc23 - 	(3*mc02*mc21+6*mc11*mc12+2*mc10*mc13+mc20*mc03+3*mc01*mc22)/mc00 + 6*(2*mc02*mc10*mc11+mc02*mc01*mc20+2*mc01*mc112)/mc002;


b3=[ mc505 mc550 mc514 mc541 mc523 mc532 ];

pq=5;
b3=[ mc505 mc550 mc514 mc541 mc523 mc532 ]/((vel(1)/2)^(pq));


% 7 8 9 10 11 12
if r>6
	mc60=mc(7,1);
	mc06=mc(1,7);
	mc51=mc(6,2);
	mc15=mc(2,6);
	mc43=mc(5,4);
	mc34=mc(4,5);
	mc33=mc(4,4);
	mc70=mc(8,1);
	mc07=mc(1,8);
	mc61=mc(7,2);
	mc16=mc(2,7);
	mc52=mc(6,3);
	mc25=mc(3,6);
	mc42=mc(5,3);
	mc24=mc(3,5);

	mc203=mc202*mc20;
	mc023=mc022*mc02;
	mc113=mc112*mc11;

	mc770=mc70 - 7*(3*mc20*mc50+5*mc40*mc30+mc10*mc60)/mc00 + 210*(mc202*mc30+mc20*mc10*mc40)/mc002 - (630*mc203*mc10)/mc003;


	mc707=mc07 - 7*(3*mc02*mc05+5*mc04*mc03+mc01*mc06)/mc00 + 210*(mc022*mc03+mc02*mc01*mc04)/mc002 - (630*mc023*mc01)/mc003;


	mc761=mc61 - (15*mc20*mc41+15*mc40*mc21+mc01*mc60+6*mc11*mc50+20*mc31*mc30+6*mc10*mc51)/mc00+30*(3*mc202*mc21+mc20*mc01*mc40+4*mc20*mc11*mc30+2*mc11*mc40*mc10+4*mc20*mc31*mc10)/mc002	- 90*(mc203*mc01+6*mc11*mc10*mc202)/mc003;

	mc716=mc16 - (15*mc02*mc14+15*mc04*mc12+mc10*mc06+6*mc11*mc05+20*mc13*mc03+6*mc01*mc15)/mc00+30*(3*mc022*mc12+mc02*mc10*mc04+4*mc02*mc11*mc03+2*mc11*mc04*mc01+4*mc02*mc13*mc01)/mc002 - 90*(mc023*mc10+6*mc11*mc01*mc022)/mc003;


	mc752= mc52 - 	(10*mc20*mc32+5*mc40*mc12+10*mc11*mc41+20*mc31*mc21+2*mc51*mc01+mc02*mc50+10*mc22*mc30+5*mc10*mc42)/mc00+10*(3*mc12*mc202+12*mc11*mc21*mc20+4*mc20*mc01*mc31+2*mc20*mc02*mc30+6*mc20*mc10*mc22+2*mc40*mc01*mc11+mc40*mc10*mc02+4*mc112*mc30+8*mc31*mc10*mc11)/mc002-90*(2*mc202*mc01*mc11+mc202*mc10*mc02+4*mc20*mc112*mc10)/mc003;

	mc725= mc25 - 	(10*mc02*mc23+5*mc04*mc21+10*mc11*mc14+20*mc13*mc12+2*mc15*mc10+mc20*mc05+10*mc22*mc03+5*mc01*mc24)/mc00+10*(3*mc21*mc022+12*mc11*mc12*mc02+4*mc02*mc10*mc13+2*mc02*mc20*mc03+6*mc02*mc01*mc22+2*mc04*mc10*mc11+mc04*mc01*mc20+4*mc112*mc03+8*mc13*mc01*mc11)/mc002-90*(2*mc022*mc10*mc11+mc022*mc01*mc20+4*mc02*mc112*mc01)/mc003;


	mc743= mc43 - 	(6*mc20*mc23+mc40*mc03+12*mc11*mc32+12*mc31*mc12+3*mc02*mc41+18*mc22*mc21+3*mc42*mc01+4*mc13*mc30+4*mc10*mc33)/mc00+6*(6*mc20*mc02*mc21+12*mc20*mc11*mc12+4*mc20*mc10*mc13+mc202*mc03+6*mc20*mc01*mc22+mc40*mc01*mc02+12*mc112*mc21+8*mc11*mc01*mc31+4*mc11*mc02*mc30+12*mc11*mc10*mc22+4*mc31*mc10*mc02)/mc002-18*(12*mc20*mc02*mc10*mc11+3*mc202*mc02*mc01+12*mc20*mc112*mc01+8*mc10*mc113)/mc003;

	mc734= mc34 - 	(6*mc02*mc32+mc04*mc30+12*mc11*mc23+12*mc13*mc21+3*mc20*mc14+18*mc22*mc12+3*mc24*mc10+4*mc31*mc03+4*mc01*mc33)/mc00+6*(6*mc02*mc20*mc12+12*mc02*mc11*mc21+4*mc02*mc01*mc31+mc022*mc30+6*mc02*mc10*mc22+mc04*mc10*mc20+12*mc112*mc12+8*mc11*mc10*mc13+4*mc11*mc20*mc03+12*mc11*mc01*mc22+4*mc13*mc01*mc20)/mc002-18*(12*mc02*mc20*mc01*mc11+3*mc022*mc20*mc10+12*mc02*mc112*mc10+8*mc01*mc113)/mc003;


b4=[mc707 mc770 mc716 mc761 mc725 mc752  mc734 mc743  ];
% 13 14 15 16 17 18 19 20
pq = 7;
b4=[mc707 mc770 mc716 mc761 mc725 mc752  mc734 mc743  ]/((vel(1)/2)^(pq));


end;
end;
end;

b=[ b1 b2 b3 b4 mc00];