function [A,U,rmax,tx,ty,x,y]= am_U_pseudonormalized(img,r,center,radius,coef) 
% [A,U,rmax,tx,ty] = am_U_pseudonormalized(img,r,center,radius,coef) computes Appell moments of the image with normalization preserving rotation invariance.
% img(n1,n2) - image matrix
% U - Appell moments
% r - 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 (implicitly 1).
% It should be set so all objects from a database were mapped into
% the unit disk.
% rmax is radius mapped onto the unit circle.
% tx,ty are coordinates mapped onto the center of the unit disk.

if nargin < 3
    center = 1;       % centroid is not stable
end
if nargin < 4
    radius = 2;
end
if nargin < 5
    coef = 1;
end

[x,y,v,rmax,tx,ty] = transform_to_unit_circle(img,center,radius,coef);

dxy = 1/rmax^2;

sx = size(x);
U{1,1} = ones(size(x))/sqrt(pi); % Appell polynomials
A = zeros(r+1);         % Appell moments
for m = 0:r/2
    if m > 0
        U{m+1,m+1} = x.*(sqrt(gamma(m+.5)/gamma(m+1))*(3*m-1)*(2*m)^(-.75)*U{m,m+1} + ...
            (m-1)*(8*m)^(-.25)*(2*m-1)^(-.75)*y.*U{m,m}); 
        if m > 1
            U{m+1,m+1} = U{m+1,m+1} + (m-1)*((y.^2*(3*m-1)-2*m+1)*m^(-1.25)*(4*m-2)^(-.75).*U{m-1,m+1}-...
                sqrt(gamma(m-.5)/gamma(m+1))*(m*(m-1)/4)^(.25)*(2*m-1)^(-.75)*y.*U{m-1,m}); 
        end
    end
    A(m+1,m+1) = v'*U{m+1,m+1}*dxy;
    for n = m:r-m-1
        [Umn_1, Um_1n_1, Um_1n, Unm_1, Un_1m, Un_1m_1] = deal(zeros(sx)); % initialize all matrices at the same time
        if m > 0
            Um_1n = U{m,n+1};
            Unm_1 = U{n+1,m};
        end
        Umn = U{m+1,n+1};
        Unm = U{n+1,m+1};
        if n > 0
            Umn_1 = U{m+1,n};
            Un_1m = U{n,m+1};
        end
        if m*n>0
            Um_1n_1 = U{m,n};
            Un_1m_1 = U{n,m};
        end
        C1 = sqrt(2)*n*(m+n+1)^(-1.25);
        C2 = m*n*(max(0,(m+n-1)))^(.25)*((m+n+1)*(max((m+n),1)))^(-.75);
        if m+n >0
            C1 = C1*(m+n)^(-.75);
            C2 = sqrt(gamma((m+n)/2)/gamma((m+n+3)/2))*C2;
        end
        C = sqrt(gamma((m+n+2)/2)/gamma((m+n+3)/2))*(m+2*n+1)*(m+n+1)^(-.75);
        U{m+1,n+2} = C*y.*Umn + C1*(m*x.*y.*Um_1n + ((x.^2-1)*n+m*(2*x.^2-1)).*Umn_1) - C2*x.*Um_1n_1;
        U{n+2,m+1} = C*x.*Unm + C1*(m*x.*y.*Unm_1 + ((y.^2-1)*n+m*(2*y.^2-1)).*Un_1m) - C2*y.*Un_1m_1;
        A(m+1,n+2) = v'*U{m+1,n+2}*dxy;
        A(n+2,m+1) = v'*U{n+2,m+1}*dxy;
    end
end
end