function [A,V,rmax,tx,ty,x,y]= am_V_pseudonormalized(img,r,center,radius,coef)
% [A,V,rmax,tx,ty] = am_pseudonormalized(img,r,center,radius,coef) computes Appell moments of the image with normalization preserving rotation invariance.
% img(n1,n2) - image matrix
% V - 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);
V{1,1} = ones(size(x))/sqrt(pi); % Appell polynomials
A = zeros(r+1);         % Appell moments
for m = 0:r/2
    if m > 0
        V{m+1,m+1} = sqrt(gamma(m+.5)/gamma(m+1))*(2*m+1)*(2*m)^(-.75)*(x.*V{m,m+1} + y.*V{m+1,m}); 
        if m > 1
            V{m+1,m+1} = V{m+1,m+1} - (m-1)*(2*m+1)*(8*m)^(-.25)*(2*m-1)^(-1.75)*(V{m-1,m+1} + V{m+1,m-1}); 
        end
    end
    A(m+1,m+1) = v'*V{m+1,m+1}*dxy;
    for n = m:r-m-1
        [Vm_2n1, Vm_1n1, Vmn_1, Vn1m_1, Vn_1m, Vn1m_2] = deal(zeros(sx)); % initialize all matrices at the same time
        if m > 0
            Vm_1n1 = V{m,n+2};
            Vn1m_1 = V{n+2,m};
            if m > 1
                Vm_2n1 = V{m-1,n+2};
                Vn1m_2 = V{n+2,m-1};
            end
        end
        Vmn = V{m+1,n+1};
        Vnm = V{n+1,m+1};
        if n > 0
            Vmn_1 = V{m+1,n};
            Vn_1m = V{n,m+1};
        end
        C1 = 2*sqrt(gamma((m+n+2)/2)/gamma((m+n+3)/2))*(m+n+1)^(-7/4);
        C2 = sqrt(2)*(max((m+n),1))^(-7/4)*(m+n+1)^(-5/4);
        V{m+1,n+2} = (m+n+2)*(C1*(x.*Vm_1n1*m + y.*Vmn*(n+1)) -... 
            C2*(Vmn_1*n*(n+1) + Vm_2n1*m*(m-1)));
        V{n+2,m+1} = (m+n+2)*(C1*(x.*Vnm*(n+1) + y.*Vn1m_1*m) -... 
            C2*(Vn_1m*n*(n+1) + Vn1m_2*m*(m-1)));
        A(m+1,n+2) = v'*V{m+1,n+2}*dxy;
        A(n+2,m+1) = v'*V{n+2,m+1}*dxy;
    end
end
end