% This program was used to generate permutation test results
% for the paper (http://www.mindboggle.info/papers/):
%
% "Evaluation of 14 nonlinear deformation algorithms 
%  applied to human brain MRI registration" 
%  NeuroImage (2009), doi:10.1016/j.neuroimage.2008.12.037
%  by Arno Klein, Jesper Andersson, Babak A. Ardekani, John Ashburner, 
%     Brian Avants, Ming-Chang Chiang, Gary E. Christensen, D. Louis Collins,
%     James Gee, Pierre Hellier, Joo Hyun Song, Mark Jenkinson, Claude Lepage,
%     Daniel Rueckert, Paul Thompson, Tom Vercauteren, Roger P. Woods, 
%     J. John Mann, and Ramin V. Parsey.
%
% (c) 2008, @rno klein
%--------------------------------------------------------------------------

nbrains = 40;

init             = 1;
run_permutations = 1;
  nruns          = 1e4; % number of times a random set of independent pairs are selected
compare_methods  = 1;
  thresh1 = 0.99;
  thresh2 = 0.0;
%s=sprintf('load output_tables/table%d.txt; table=table%d;',nbrains,nbrains); eval(s);
s=sprintf('load output_tables/table%d_avgTO.txt; table=table%d_avgTO;',nbrains,nbrains); eval(s);
All_brain_pairs  = table(:,1:2);
All_brain_values = table(:,3:end);
nmethods         = size(All_brain_values,2);
method_names = {'FLIRT','AIR','ANIMAL','ART','D.Demons','FNIRT','IRTK','JRD-fluid','ROMEO','SICLE','SyN','SPM_N*','SPM_N','SPM_US','SPM_D*'};%,'SPM_D'};

%------------------------------

nbrain_pairs     = nbrains/2;
method_pairs     = nchoosek(1:nmethods,2);
nmethod_pairs    = size(method_pairs,1);

% number of times columns are randomly swapped
if nbrains > 12, 
     nmixall = 0;
     nmixs   = 1000;
     nskip   = 1;
     
else nskip = 1;
     nmixall = 1;
end

if init==1
 if nbrains<=12
   p1 = nchoosek([1:nbrains],nbrain_pairs);
   p2 = unique(p1>nbrain_pairs,'rows');
   p3 = perms(1:nbrain_pairs);
   for ip = 1:size(p3,1)
       p0 = unique(p1(:,p3(ip,:))>nbrain_pairs,'rows');
       p2 = [p2; p0];
   end
   binary_perms = unique(p2,'rows')';
 else
   binary_perms = (round(rand(nmixs,nbrain_pairs))>0.5)';
 end
 s=sprintf('save output_tables/binaryperms%d_avgTO.mat binary_perms;',nbrains); eval(s);
end

if run_permutations==1
    
    s=sprintf('load output_tables/binaryperms%d_avgTO.mat;',nbrains); eval(s);
    if nmixall==1
       nmixs = size(binary_perms,2)-1;
    end
    
    ZD  = zeros(nmixs,1);
    pvalues = zeros(nruns,nmethod_pairs);
    
    for irun = 1:nruns

        PERM_brain_pairs = randperm(nbrains);
        PERM_brain_pairs = reshape(PERM_brain_pairs,nbrain_pairs,2);
        V                = zeros(nbrain_pairs,nmethods);

        for ibrain_pair=1:nbrain_pairs
            I=find(All_brain_pairs(:,1)==PERM_brain_pairs(ibrain_pair,1) & All_brain_pairs(:,2)==PERM_brain_pairs(ibrain_pair,2));
            if size(I,1)>0
                V(ibrain_pair,:) = All_brain_values(I(1),:);
            end
        end

        for imethod_pair = 1:nmethod_pairs
          
          if sum(V(:,method_pairs(imethod_pair,:)))>0  
              VP = V(:,method_pairs(imethod_pair,:));
              D0 = sum(VP(:,1) - VP(:,2))/nbrain_pairs;
              D  = ZD;
              for imix=2:nskip:nmixs
                I = binary_perms(:,imix);
                M = [VP(:,1).*I + VP(:,2).*(1-I) ...
                     VP(:,2).*I + VP(:,1).*(1-I)];
                D(imix) = sum(M(:,1) - M(:,2))/nbrain_pairs;
              end

              I  = find(abs(D)>abs(D0));
              pvalues(irun,imethod_pair) = (1 + size(I,1)) / (nmixs+1);
          end
        end
        
        irun

    end
    
    %pvalues
    num_pvalues_05  = sum(pvalues<0.05&pvalues>0,1);
    perc_pvalues_05 = num_pvalues_05/nruns;
    
    s=sprintf('save output_tables/permpvalues%d_avgTO.mat *pvalues* nruns nmixs method_pairs nmethod_pairs method_names;',nbrains); eval(s);  

end

if compare_methods==1

   MEANS = mean(All_brain_values,1);
   s=sprintf('load output_tables/permpvalues%d_avgTO.mat;',nbrains); eval(s);  
   M1 = zeros(nmethods,nmethods);
   M2 = M1;
   for i=1:size(method_pairs,1)
       if perc_pvalues_05(i) >= thresh1
          if MEANS(method_pairs(i,1))>MEANS(method_pairs(i,2))
             M1(method_pairs(i,1),method_pairs(i,2)) =  perc_pvalues_05(i);
          else
             M1(method_pairs(i,1),method_pairs(i,2)) = -perc_pvalues_05(i);
          end
       end
       if perc_pvalues_05(i) >= thresh2
           
          if MEANS(method_pairs(i,1))>MEANS(method_pairs(i,2))
             M2(method_pairs(i,1),method_pairs(i,2)) =  perc_pvalues_05(i);
          else
             M2(method_pairs(i,1),method_pairs(i,2)) = -perc_pvalues_05(i);
          end
       end
   end
   
   M1 = M1 - M1';
   M2 = M2 - M2';
  
   subplot(1,2,1)
   imagesc(M1); axis square; colorbar
   set(gca,'YTick',[1:nmethods],'YTickLabel',method_names,'FontSize',10,'FontWeight','Bold')
   set(gca,'XTick',[1:nmethods],'XTickLabel',method_names,'FontSize',10,'FontWeight','Bold')
   %xticklabel_rotate([1:nmethods],60,char(method_names),'interpreter','none')
   s=sprintf('title(''Threshold: %1.2f'');',thresh1); eval(s);
   
   subplot(1,2,2)
   imagesc(M2); axis square; colorbar
   set(gca,'YTick',[1:nmethods],'YTickLabel',method_names,'FontSize',10,'FontWeight','Bold')
   set(gca,'XTick',[1:nmethods],'XTickLabel',method_names,'FontSize',10,'FontWeight','Bold')
   %xticklabel_rotate([1:nmethods],60,char(method_names),'interpreter','none')
   s=sprintf('title(''Threshold: %1.2f'');',thresh2); eval(s);

   mean_perc_pvalues = mean(M1,2);
   std_perc_pvalues  = std(M1,[],2);
   [sortmeans,I] = sort(mean_perc_pvalues,'descend');
   disp('Threshold 1');
   for i=1:nmethods
     s=sprintf('disp(''%1.2f %1.2f %s'');',sortmeans(i),std_perc_pvalues(I(i)),char(method_names(I(i)))); eval(s);
   end

   mean_perc_pvalues = mean(M2,2);
   std_perc_pvalues  = std(M2,[],2);
   [sortmeans,I] = sort(mean_perc_pvalues,'descend');
   disp('Threshold 2');
   for i=1:nmethods
     s=sprintf('disp(''%1.2f %1.2f %s'');',sortmeans(i),std_perc_pvalues(I(i)),char(method_names(I(i)))); eval(s);
   end
   
   disp('means means+SDmax means+SDmax*2 method');
   SDmax = std_perc_pvalues(I(1));
   for i=1:nmethods
     s=sprintf('disp(''%1.2f %1.2f %1.2f %s'');',...
       sortmeans(i),sortmeans(i)+SDmax,sortmeans(i)+2*SDmax,char(method_names(I(i)))); eval(s);
   end

end