function [ phi, ift ] = RotInvsSym(CM,r,t,thres)
%normalizace k rotaci
% Hledá se normalizační moment s dostatečnou absolutní hodnotou schopný
% normalizovat i symetrické obrázky.
%
% CM matice momentů
% r maximální řád momentů
% t typ normalizace k rotaci
% t=1; odečte se příslušný násobek fáze normalizačního momentu
% t=2; vynásobí se příslušnou mocninou normalizačního momentu
% thres práh absolutní hodnoty normalizačního momentu
% output:
% phi is matrix with the invariants, 
% ift vecotor with information: ift=[p0,q0,abs(CM(p0+1,q0+1))]

if nargin<2
    r=min(size(CM))-1;
end
if nargin<3
    t=2;
end
if nargin<4
    thres=1e-3;
end

flag=1;
p0=-1;
q0=-1;
id=0;   %index difference
while id<r && flag==1
    id=id+1;
    p=0;    %first index
    if id==1
        p=1;
    end
    while p<=floor((r-id)/2) && flag==1
        q=p+id;
        if abs(CM(p+1,q+1))>thres
            flag=0;
            p0=p;
            q0=q;
        end
        p=p+1;
    end
end

[X,Y] = meshgrid(0:r);
if flag==0
%     disp(['normalizační moment A_',int2str(p0),',',int2str(q0),' absolutní hodnota: ',num2str(abs(CM(p0+1,q0+1)))])
    if t==2
        A = CM(p0+1,q0+1).^((Y-X)/id);
    elseif t==1
        Theta = angle(CM)*(id-1)+angle(CM(p0+1,q0+1))*(Y-X);
%         Theta = angle(CM(p0+1,q0+1))*((Y-X)/id);
%         Theta = normangle(angle(CM(p0+1,q0+1))*(Y-X))/id;        
        A = exp(1i*Theta);
    end
else
%     disp('žádný nadprahový normalizační moment se nenašel')
    A = ones(r+1);
end
N = Y>=X;
% phi = CM.*A;
phi = CM.*A.*N;
% phi = CM.*N;
% if flag==0
%     disp([abs(CM(5,1)),angle(CM(5,1))*180/pi,abs(CM(p0+1,q0+1)),angle(CM(p0+1,q0+1))*180/pi,...
%           p0,q0,abs(phi(5,1)),angle(phi(5,1))*180/pi,(Y(5,1)-X(5,1))/id])
% else
%     disp([abs(CM(5,1)),angle(CM(5,1))*180/pi,1,0,-1,-1,abs(phi(5,1)),angle(phi(5,1))*180/pi])
% end
if flag==0
    ift=[p0,q0,abs(CM(p0+1,q0+1))];
else
    ift=[p0,q0,0];
end
end