
%% set up parameters

% set run times
run=100;  
generate_run=cell(1,run);

recP=[]; 
% set sythesized image size
xx=500; yy=500;  

% set synthetic transformations P
 for i=1:run+500  
  a=500; h=-1000;  
  P=space_generateP(a,h);    
  generate_run{i}.P=P;  
 end

% set feature point numbers std_n
   n=8000;           
   std_n=200;
   
% set outliers percentage
   outliers_pct=[0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8];
   noise_degree=1;    


%% run each random trial everytime
for t=1:run 
    t
    P=generate_run{t}.P;     
%% 
    
 % generate and adjust matched feature point sets
     X=[];Y=[];ii=0; 
     while length(X)<(std_n+200) || length(Y)<(std_n+200) 
         ii=ii+1;
         if ii>=2;P=generate_run{run+floor(500*rand(1))}.P;end
         X=[a*2*(rand(1,n)-0.5);a*2*(rand(1,n)-0.5);ones(1,n)] ;

%%
         Y=P*X;  Y(1,:)=Y(1,:)./Y(3,:); Y(2,:)=Y(2,:)./Y(3,:); Y(3,:)=1;

         x_runout=union(find(Y(1,:)>xx),find(Y(1,:)<-xx));   
         y_runout=union(find(Y(2,:)>yy),find(Y(2,:)<-yy));   
         Yrunout=union(x_runout,y_runout);
         Y=Y(:,setdiff(1:n,Yrunout));
         X=X(:,setdiff(1:n,Yrunout));
     end

     kk=size(Y,2);S=randperm(kk,min(kk,std_n));
     X=X(:,S);
     Y=Y(:,S);
     Xini=X;  Yini=Y;
     rec_compare3=cell(size(outliers_pct,2),size(noise_degree,2));
 % add outliers
        for i=1:size(outliers_pct,2)
          for j=1:size(noise_degree,2)
              X=Xini;  Y=Yini; 
              std_n=length(X);
              Y(1:2,:)=Y(1:2,:)+normrnd(0,1,2,size(Y,2));
              rand_pos=randperm(std_n,floor(outliers_pct(i)*std_n)) ;
              Y(1,rand_pos)=2*a*(rand(1,size(rand_pos,2))-0.5);
              Y(2,rand_pos)=2*a*(rand(1,size(rand_pos,2))-0.5);
 % get index of synthesized right inliers
              K=P*X; K(1,:)=K(1,:)./K(3,:); K(2,:)=K(2,:)./K(3,:); K(3,:)=1;
              V=Y-K;                                                                       
              D=V(1,:).^2+V(2,:).^2;                                             
              real_in=sort(find(D<=25));
              
 % compare and record
               tic;
                  recoverin_AHC=AHC(X,Y,5);           
               timecost_AHC=toc;

                 rec_compare3{i,j}.timecost.AHC=timecost_AHC;
                 rec_compare3{i,j}.P.AHC=length(intersect(real_in,recoverin_AHC))/length(recoverin_AHC);
                 rec_compare3{i,j}.R.AHC=length(intersect(real_in,recoverin_AHC))/length(real_in);
                 rec_compare3{i,j}.F.AHC=2*rec_compare3{i,j}.P.AHC*rec_compare3{i,j}.R.AHC/...
                 (rec_compare3{i,j}.P.AHC+rec_compare3{i,j}.R.AHC);
        
               if size(intersect(real_in,recoverin_AHC),2)==0 || isnan(rec_compare3{i,j}.F.AHC)
               rec_compare3{i,j}.P.AHC=0;  rec_compare3{i,j}.R.AHC=0; rec_compare3{i,j}.F.AHC=0; 
               end

          end
        end 
 
       generate_run{t}.recPRFT=rec_compare3;

end 


%% calculate the average of time, precision score, recall score and F-score
AHC_avgT=zeros(size(outliers_pct,2),size(noise_degree,2));
AHC_avgP=zeros(size(outliers_pct,2),size(noise_degree,2));
AHC_avgR=zeros(size(outliers_pct,2),size(noise_degree,2));
AHC_avgF=zeros(size(outliers_pct,2),size(noise_degree,2));


for t=1:run
             for i=1:size(outliers_pct,2)
                for j=1:size(noise_degree,2)
                    AHC_avgT(i,j)=AHC_avgT(i,j)+generate_run{t}.recPRFT{i,j}.timecost.AHC;
                    AHC_avgP(i,j)=AHC_avgP(i,j)+generate_run{t}.recPRFT{i,j}.P.AHC;
                    AHC_avgR(i,j)=AHC_avgR(i,j)+generate_run{t}.recPRFT{i,j}.R.AHC;
                    AHC_avgF(i,j)=AHC_avgF(i,j)+generate_run{t}.recPRFT{i,j}.F.AHC;    
                      
                end
             end 
             
end 

AHC_avgT=AHC_avgT/run;
AHC_avgP=AHC_avgP/run;
AHC_avgR=AHC_avgR/run;
AHC_avgF=AHC_avgF/run;

               

%% Figure the result               
    title('Comparison at Different Levels of Outliers (1000 pairs)');

    ylabel('Average F-score');xlabel('Percentage of Outliers under Projective transformation');
    axis([1,8,0.1,1]);
    set(gca,'xtick',1:1:8);
    set(gca,'xticklabel',{'10','20','30','40','50','60','70','80'});


    grid on;
    i=1:8;
    hold on;
    plot(i,AHC_avgF(i,1),'Marker','o','Color','r');
    hold on;
    legend('AHC');
    hold off;