function [invts,invmat,ind,normom,idcoef]=chebyshev_fourier(a,r,kind,thres,center,radius,coef)
% [invts,invmat,ind,normom]=chebushev_fourier(a,r,kind,thres) calculates 
% values of rotation Chebushev-Fourier invariants of the image a (no impl. value) 
% up to the order r (impl. 3), 
% kind is kind of the Chebyshev polynomial, 1 or 2,
% thres is threshold of non-zero moments, impl. 1e-3
% 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.

% invts  - values of the invariants
% invmat - values of the invariants in matrix form,
%            A_{r,ell}=invmat(r+1,(r+ell)/2+1);
% ind    - vector of indicators, ind(1,i)=i is the number of i-th invariant,
%          ind(2,i) is order of i-th invariant, ind(3,i) is repetition and
%          ind(4,i)=0 for real part and ind(4,i)=1 for imaginary part.
% normom - vector of three numbers describing the moment normalizing rotation:
%          order, repetition and normalizing angle.

if nargin<2
    r=3;
end
if nargin<3
    kind=1;
end
if nargin<4
    thres=1e-3;
end

[x1,y1]=find(a);
szx=size(x1);
szy=size(y1);
if(szx(1)>0 && szy(1)>0)
	mn1=min(x1);
	mx1=max(x1);
	mn2=min(y1);
	mx2=max(y1);
    a=a(mn1:mx1,mn2:mx2);
end

a=a/max(a(:));  %normalization of moments to contrast
[am,~,~,~,rma]=chfm(a,r,kind,center,radius,coef);   %computation of Chebyshev-Fourier moments
idcoef=rma(2)/rma(3);
mt=0;
nt=0;
theta=0;
flag=0;         %normalization of moments to rotation
for k1=1:r
    for k2=k1:2:r
        if k2>1
            if abs(am(k2+1,(k2+k1)/2+1))>thres
                theta=angle(am(k2+1,(k2+k1)/2+1))/k1;   %normalizing moment found
                flag=1;
                mt=k2;
                nt=k1;
            end
        end
        if flag
            break
        end
    end
    if flag
        break
    end
end
if flag
    for k1=2:r
        for k2=-k1:2:k1
            am(k1+1,(k1+k2)/2+1)=am(k1+1,(k1+k2)/2+1)*exp(-1i*k2*theta);
        end
    end
end
%end of normalization of moments to rotation

% invariants
pinv=1;
for k1=2:r
    for k2=k1:-2:0
        invts(pinv)=real(am(k1+1,(k1+k2)/2+1));
        cmp(pinv,:)=[pinv,k1,k2,0];
        pinv=pinv+1;
        if k2~=0 && (k1~=3 || k2~=1)
            invts(pinv)=imag(am(k1+1,(k1+k2)/2+1));
            cmp(pinv,:)=[pinv,k1,k2,1];
            pinv=pinv+1;
        end
    end
end
invmat=am;
ind=cmp';
normom=[mt,nt,theta];