%experiment with database ALOI.

clear invs c

tic

normmom=4;  %normmom=0; similarity invariants from complex moments, 
            %normmom=1; Zernike moments, 
            %normmom=2; Gauss-Hermite moments, 
            %normmom=3; Chebyshev-Fourier moments, 
            %normmom=4; Appell moments, 

rr = 5;     % maximum order
nr = 36;    % the number of rotations
nc = 1000;  % the number of classes (ALOI max. 1000)
p = 'U';    % p='U' for Appell U, p='V' for Appell V, 
tr = 1;     % tr type of normalization to rotation of Appell moments
            % tr=1; subtraction of phase c_12
            % tr=2; multiplication by c_12
switch p
    case 'U'
        th=0.0001; %threshold of Appell moments
    case 'V'
        th=0.01;   %threshold of Appell moments
end
center = 2; % center=1 center, center=2 centroid of the image is mapped onto the center of the unit disk.
radius = 2; % 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 = 1;   % 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.

dir= 'F:\images\ALOI\aloi_grey_view\';   %directory with dataset
drm= 'F:\images\ALOI\aloi_mask\';        %directory with masks
            
sigma=0.3;  %standard deviation parameter of the GH polynomials
kind=1;     %kind of Chebyshev polynomials, 1 or 2
thres=1e-3; %threshold of Zernike moments
typec=0;    %type is method of computation of Zernike moments: 
            %if typec=0, moments are computed directly else they are converted from geometric ones
typei=0;    %if typei=0, rotation normalization is used else phase cancellation by multiplication is used
            
snrk=1000;    %PSNR of the Gaussian noise in decibells

seed=0;
rng('default'); %random sequence is the same in every time

pf=(rr+1)*(rr+2)/2-4;
if normmom==4
    indices = zeros(pf,3);
    np=0;
    for i=0:rr
        for j=0:i/2
            m=i-j;
            n=j;
            np=np+1;
            indices(np,:)=[m,n,0];
            if n<m 
                np=np+1;
                indices(np,:)=[m,n,1];
            end
        end
    end
    pf=np;
end

h = waitbar(0,'1','Name','ALOI recognition','CreateCancelBtn','setappdata(gcbf,''canceling'',1)');
setappdata(h,'canceling',0);

invs=zeros(nr*nc,pf);
pat=zeros(nr*nc,pf);
qat=zeros(nr*nc,pf);
pa=zeros(1,pf);
qa=zeros(1,pf);
ica=zeros(1,nc);
imga=zeros(5760,7680);  % image with examples of ALOI dataset
imgra=zeros(3000,3000);  % image with rotations of an example from ALOI dataset
for i=1:nc
    if getappdata(h,'canceling')
        break
    end
    waitbar((i-1)/nc,h,sprintf('features of class %d from %d',i,nc))
    
    if i==60
        t=0;
    end
    img=double(imread([dir,int2str(i),'\',int2str(i),'_r0.png']))/255;
    msk=double(imread([drm,int2str(i),'\',int2str(i),'_c1.png']));
    img=img.*msk;
    if max(img(:))==0
        disp(['zero image ',int2str(i)])
    end
    if mod(i,11)==1
        ia=floor((i-1)/11/10);
        ja=mod((i-1)/11,10);
        imga(ia*size(img,1)+1:(ia+1)*size(img,1),ja*size(img,2)+1:(ja+1)*size(img,2))=img;
    end
    for t=0:nr-1
        if t>0
            imgr=imrotate(img,round(t*360/nr),'nearest','loose');
            mskr=imrotate(msk,round(t*360/nr),'nearest','loose');
        else
            imgr=img;
            mskr=msk;
        end
        if snrk>=0 && snrk<1000
            stds=std(imgr(:));
            nstd=10^(-snrk/20)*stds;
            imgr=imgr+nstd*mskr.*randn(size(imgr));
        end
        if i==155
            [yc,xc]=find(imgr);
            mn1=min(xc);
            mx1=max(xc);
            mn2=min(yc);
            mx2=max(yc);
            imgrb=imgr(mn2:mx2,mn1:mx1);
            if t==0
                [nr1,nr2]=size(imgrb);
                nrc=max(nr1,nr2);
                nrc=round(nrc*sqrt(2));
            end
            ia=floor(t/6)*nrc+round((nrc-size(imgrb,1))/2);
            ja=mod(t,6)*nrc+round((nrc-size(imgrb,2))/2);
            imgra(ia+1:ia+size(imgrb,1),ja+1:ja+size(imgrb,2))=imgrb;
        end
        if normmom==0
            c=rotmi(imgr,rr);
        elseif normmom==1
            c=zermi(imgr,rr,thres,typec,typei,center,radius,coef);
        elseif normmom==2
            c=gauss_hermite(imgr,rr,sigma);
        elseif normmom==3
            [c,~,~,~,idcoef]=chebyshev_fourier(imgr,rr,kind,thres,center,radius,coef);
            if t==0
                ica(i)=idcoef;
            end
        elseif normmom==4
            [cm,ccm,ift]=Appell_rotation_invariants(imgr,rr,p,tr,th,center,radius,coef);
            for k=1:pf
                mf=indices(k,1);
                nf=indices(k,2);
                if indices(k,3)==0
                    c(k)=real(cm(mf+1,nf+1));
                else
                    c(k)=imag(cm(mf+1,nf+1));
                end
                pa(k)=ift(1);
                qa(k)=ift(2);
            end
            pat((i-1)*nr+t+1,:)=pa;
            qat((i-1)*nr+t+1,:)=qa;
        end
        invs((i-1)*nr+t+1,:)=c;
    end
end

eri=zeros(1,nc);
for i=1:nc
    for j=1:size(invs,2)
        mv=mean(invs((i-1)*nr+1:(i-1)*nr+nr,j));
        erc=0;
        if abs(mv)>1e-3
            for t=1:nr
                c=abs((invs((i-1)*nr+t,j)-mv)/mv);
                erc=erc+c;
            end
        end
        erc=erc/nr;
        eri(i)=eri(i)+erc;
    end
    eri(i)=eri(i)/size(invs,2);
end
err=mean(eri);
disp(['Average relative error of the invariants is ',num2str(err*100),'%'])

cha=zeros(1,nr/4);
mcca=zeros(1,nr/4);
mcci=zeros(nc,nr/4);
for k=1:nr/4
    if getappdata(h,'canceling')
        break
    end
    waitbar(4*k/nr,h,sprintf('rotation %d from %d is class representative',k,nr/4))
    ch=0;
    ctg=zeros(nc);
    for i=1:nc
        mcc=0;
        ncc=0;
        for j=1:nr-1
            j1=j;
            if j>=k
                j1=j+1;
            end
            dif=zeros(1,nc);
            c=invs((i-1)*nr+j1,:);
            for i1=1:nc
                c1=invs((i1-1)*nr+k,:);
                dif(i1)=sqrt((c-c1)*(c-c1)');
            end
            [mna,mni]=sort(dif);
            mn=mna(1);
            pmn=mni(1);
            ctg(pmn,i)=ctg(pmn,i)+1;
            if pmn~=i
                ch=ch+1;
            else
                if mna(1)>0
                    cc=mna(2)/mna(1);
                    mcc=mcc+cc;
                    ncc=ncc+1;
                end
            end
        end
        if ncc>0
            mcci(i,k)=mcc/ncc;
        end
    end
    disp(['The number of errors for rotation ',num2str((k-1)*360/nr),'° as class representative is ',int2str(ch)]);
    cha(k)=ch;
end

delete(h)

disp('confidence coefficients - mean values')
for k=1:nr/4
    for i=1:nc
        mcca(k)=mcca(k)+mcci(i,k)*ctg(i,i);
    end
    mcca(k)=mcca(k)/sum(diag(ctg));
    disp(mcca(k))
end

toc