function [A,rmax,tx,ty,rma]=chfm(img,rd,kind,center,radius,coef)
% [A]=chfm(img,rd,center,radius,coef) computes Chebyshev-Fourier moments of the image.
% img(n1,n2) - image matrix
% rd - maximum order
% kind - kind of the Chebyshev polynomial, 1 or 2.
% 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.
% rma - array with all possible rmax, rmax(i) is rmax for radius=i with
% coef=1.

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

[n1,n2]=size(img);

tx=(n2+1)/2;
ty=(n1+1)/2;
m00 =sum(img(:));
mc=max(img(:));
if center==2
	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
    %sqrt(n2/n1+n1/n2)/sqrt(2) is correction for rectangular images
    %rmax=sqrt(m00/mc)*sqrt(n2/n1+n1/n2);   
    rmax=sqrt(m00/mc);   
end

rma=zeros(1,3);
rma(1)=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]));
rma(2)=max((x.^2+y.^2).^0.5);
rma(3)=sqrt(m00/mc)*sqrt(n2/n1+n1/n2);   

rmax=rmax*coef;

%disp([rmax,tx,ty])

x=x/rmax;
y=y/rmax;
r=sqrt(x.^2+y.^2);
theta=atan2(y,x);
idx=find(abs(r)<=1);
r=r(idx);
theta=theta(idx);
v=v(idx);

norm=zeros(1,rd+1);
norm(1)=pi/(7*kind-6);
norm(2)=pi/(6*kind-4);

n2=length(r);
rm=zeros(rd+1,n2);
rm(1,:)=ones(1,n2);  % polynomial of order 0
if rd>0
    rm(2,:)=kind*(2*r'-ones(1,n2)); % polynomial of order 1: kind*x = kind*(2r-1)
end
for m=3:rd+1 %order
    rm(m,:)=(4*r'-2).*rm(m-1,:)-rm(m-2,:); % 2*(2r-1)*rm(m-1,:)-rm(m-2,:)
    norm(m)=pi/(6*kind-4);
end

weight=real((4*(r-r.^2)).^(kind-1.5));
weight(isinf(weight))=0;
weight(isnan(weight))=0;
weight(weight>10)=10;
weight(weight<-10)=-10;

for n=0:rd
    rm(n+1,:)=rm(n+1,:).*sqrt(weight'/norm(n+1));
end

for m=0:rd %order
    for n=-m:2:m  %repetition
        vmn=v'.*rm(m+1,:).*exp(-1i*n*theta');
        A(m+1,(m+n)/2+1)=sum(vmn(:));
    end
end
A=A/rmax^2;