function [M,tx,ty]=cm(l,r)

% M is a matrix of central moments up to the r-th order of the image l
% l(n1,n2) - image matrix
% tx, ty coordinates of the centroid
% The moment \mu_{pq} = M(p+1,q+1)

l=double(l);
[n1,n2]=size(l);

m00 =sum(sum(l));
w=1:n2;
v=1:n1;
if m00~=0
    tx=(sum(l*w'))/m00; 
    ty=(sum(v*l))/m00;
else
    tx=(n2+1)/2;
    ty=(n1+1)/2;
end

a=w-tx;
c=v-ty;

A=repmat(a,r+1,1).^repmat((0:r)',1,n2);
C=repmat(c,r+1,1).^repmat((0:r)',1,n1);
M=A*l'*C';

if r>0
    M(1,2)=0;
    M(2,1)=0;
end;