function [invts,pqi,aspect]=gauss_hermite(a,mxr,sigma,normcoef,typemap,mapcoef)
% [invts,pqi]=gauss_hermite(a,mxr,sigma) computes values of 
% translation and rotation invariants based on Gaussian - Hermite moments.
% a is input image,
% mxr is maximum order of invariants,
% sigma is standard deviation parameter of the GH polynomials
% invts are the invariants
% pqi: pqi(1,ni) is 1st index, pqi(2,ni) is 2nd index of ni-th invariant, 
% if pqi(3,ni)==0, the ni-th invariants is based on the real part of 
% the moment GH_{pqi(1,ni),pqi(2,ni)} else on the imaginary part.
% normcoef is normalizing coefficient.
% normcoef=1, ghm_{p0} and ghm_{0q} are normalized precisely,
% normcoef=0, ghm_{p/2,p/2} are normalized precisely. Default is 0.5.
% typemap is type of mapping.
% typemap=1, mapping according to the most distant non-zero pixel
% typemap=2, mapping according to m00
% typemap=3, mapping according to the size of the image
% typemap=4, mapping according to fixed size of images in mapcoef
% mapcoef: the mapping is multipled by mapcoef

if nargin<3
%    sigma=0.9018*(1.3218*mxr+12.734)^(-0.521);
    sigma=0.3;
end
if nargin<4
    normcoef=0.5;
end
if normcoef<0
    normcoef=0;
end
if normcoef>1
    normcoef=1;
end
if nargin<5
    typemap=1;
end
if nargin<6
    mapcoef=1;
end

% range of input values of the Gaussian - Hermite polynomial is from
%-1/sigma to 1/sigma

if mxr==2
    p0=0;
    q0=2;
else
    p0=1;
    q0=2;
end
[n1,n2]=size(a);
m00 =sum(a(:));
w=1:n2;
v=1:n1;
if m00~=0
    tx=(sum(a*w'))/m00; 
    ty=(sum(v*a))/m00;
else
    tx=(n2+1)/2;
    ty=(n1+1)/2;
end
tx=round(tx);
ty=round(ty);
[yv,xv]=find(a);
if isempty(xv) || isempty (yv)
    xv=[1;1;n2;n2];
    yv=[1;n1;n1;1];
end    
lnp=length(xv);
mxxy1=max(sqrt((xv-tx.*ones(lnp,1)).^2+(yv-ty.*ones(lnp,1)).^2));
mxxy2=sqrt(m00)/sqrt(pi/2);
aspect=mxxy1/mxxy2;
if typemap==1
    mxxy=mxxy1;
elseif typemap==2
    mxxy=mxxy2;
elseif typemap==3
    xv=[1;1;n2;n2];
    yv=[1;n1;n1;1];
    lnp=length(xv);
    mxxy=max(sqrt((xv-tx.*ones(lnp,1)).^2+(yv-ty.*ones(lnp,1)).^2));
else
    mxxy=1;
end
mxxy=ceil(mxxy*mapcoef);
mxxy1=ceil(mxxy1);
xy=(-mxxy1:mxxy1)/(sigma*mxxy);

ghvalues=zeros(mxr+1,2*mxxy1+1);
ghvalues(1,:)=exp(-xy.^2/2);
ghvalues(2,:)=2*xy.*exp(-xy.^2/2);
for n=2:mxr
    ghvalues(n+1,:)=2*xy.*ghvalues(n,:)-(2*n-2)*ghvalues(n-1,:);
end

% Gaussian - Hermite moments
ghm=zeros(mxr+1,mxr+1,3);
c=zeros(mxr+1,mxr+1,3);
xi=(1:n2)-tx+mxxy1+1;
xi(xi<1)=1;
xi(xi>2*mxxy1+1)=2*mxxy1+1;
yi=(1:n1)-ty+mxxy1+1;
yi(yi<1)=1;
yi(yi>2*mxxy1+1)=2*mxxy1+1;
for p=0:mxr
    for q=0:mxr-p
        xgh=ghvalues(p+1,xi);
        ygh=ghvalues(q+1,yi);
        ghm(p+1,q+1)=ygh*a*xgh'/mxxy^2;
%        ghm(p+1,q+1)=ghm(p+1,q+1)*(2^(p+q)*factorial(p+q)*pi^0.5)^(-0.5);
        ghm(p+1,q+1)=ghm(p+1,q+1)*(2^(p+q)*gamma(p+q+1)^normcoef*gamma((p+q)/2+1)^(2*(1-normcoef))*pi)^(-0.5);
    end
end

c=cmfromgm(mxr,ghm);                   %conversion of the central moments to complex
% rotation invariants
id=q0-p0;                               %index difference
ni=1;                                   %sequential number of the invariant
for r1=max(2,id):id:mxr
    for p=round(r1/2):r1
        q=r1-p;
        if mod(p-q,id)==0
            invts(ni)=real(c(p+1,q+1)*c(p0+1,q0+1)^((p-q)/id));
            pwi(ni)=1+(p-q)/id;
            pqi(1,ni)=p;
            pqi(2,ni)=q;
            pqi(3,ni)=0;
            ni=ni+1;
%            sprintf('Re(c_{%d%d}c_{%d%d}^%d)',p,q,p0,q0,(p-q)/id)
            if p>q && (p~=q0 || q~=p0)
%                invts(ni)=imag(c(p+1,q+1)*c(p0+1,q0+1)^((p-q)/id));
                invts(ni)=real(c(p+1,q+1)*c(q+1,p+1));
%                pwi(ni)=1+(p-q)/id;
                pwi(ni)=2;
                pqi(1,ni)=p;
                pqi(2,ni)=q;
                pqi(3,ni)=1;
                ni=ni+1;
%                sprintf('Im(c_{%d%d}c_{%d%d}^%d)',p,q,p0,q0,(p-q)/id)
            end
        end
    end
end
invts=sign(invts).*abs(invts).^(1./pwi);    %magnitude normalization to degree