function [A,rmax,tx,ty]=zm(img,rd,center,radius,coef)
% [A]=zm(img,rd,center,radius,coef) computes Zernike moments of the image.
% img(n1,n2) - image matrix
% rd - maximum order
% If center=1 center of the image is mapped onto the center of the unit disk.
% If center=2 centroid of the image is mapped onto the center of the unit disk.
% If radius=1, the most distant corner is mapped onto the unit circle.
% If radius=2, the most distant nonzero pixel is mapped onto the unit circle.
% If radius=3, the sqrt(m00), is mapped onto the unit circle.
% coef - the radius mapped onto the unit circle is multiplied by coef (implicitely 1). 
% It should be set so all objects from a database were mapped into 
% the unit disk.
% The moment A_{mn} = A(m+1,(m+n)/2+1).
% rmax is radius mapped onto the unit circle.
% tx,ty are coordinates mapped onto the center of the unit disk.

if nargin<3
    center=2;
end
if nargin<4
    radius=3;
end
if nargin<5
    coef=1;
end

img=double(img);
[n1,n2,n3]=size(img);
if n3>1
    img=sum(img,3)/n3;
end


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

rmax=sqrt(max([(1-tx)^2+(1-ty)^2,(n2-tx)^2+(1-ty)^2,(n2-tx)^2+(n1-ty)^2,(1-tx)^2+(n1-ty)^2]));
A=zeros(rd+1,rd+1);
[y,x,v]=find(img);
if isempty(v)
    return
end
x=x-tx;
y=y-ty;
if radius==2
    rmax=max((x.^2+y.^2).^0.5);
elseif radius==3
    mc=max(img(:));
    %sqrt(n2/n1+n1/n2)/sqrt(2) is correction for rectangular images
    rmax=sqrt(m00/mc)*sqrt(n2/n1+n1/n2)/sqrt(2);   
end
rmax=rmax*coef;

%disp([rmax,tx,ty])

x=x/rmax;
y=y/rmax;
r=sqrt(x.^2+y.^2);
theta=atan2(y,x);

%Kintner method
for n=-rd:rd  %repetition
    an=abs(n);
    rmn0=v.*r.^an;
    vmn=rmn0.*exp(-1i*n*theta);
    A(an+1,(an+n)/2+1)=(an+1)/pi*sum(vmn(:));
    if(rd-an>=2)
        rmn2=v.*((an+2)*r.^(an+2)-(an+1)*r.^an);
        vmn=rmn2.*exp(-1i*n*theta);
        A(an+3,(an+2+n)/2+1)=(an+3)/pi*sum(vmn(:));
    end
    for m=an+4:2:rd %order
        k1=(m+n)*(m-n)*(m-2)/2;
        k2=2*m*(m-1)*(m-2);
        k3=-n^2*(m-1)-m*(m-1)*(m-2);
        k4=-m*(m+n-2)*(m-n-2)/2;
        rmn4=((k2*r.^2+k3*ones(size(v))).*rmn2+k4*rmn0)/k1;
        vmn=rmn4.*exp(-1i*n*theta);
        A(m+1,(m+n)/2+1)=(m+1)/pi*sum(vmn(:));
        rmn0=rmn2;
        rmn2=rmn4;
    end
end
A=A/rmax^2;