% This program was used to generate some of the results and figures
% 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.
%
% This program compares source-target brain registration results 
% across different registration methods. 
% The results were generated from evaluate_warps_main.m and 
% evaluate_distances_main.m and compares manual labels for the source brain 
% registered to manual labels of the target brain, using different measures:
%
% * Target, union, or mean overlap for volumes or surfaces 
% * Volume similarity 
% * Average distance error between label boundaries 
%
% These results may be computed for each brain region, 
% or for the entire brain, averaged across regions.
%
% Results take the form of plots (box plots, line plots, color tables)
% or tables of values (corresponding to the color tables, 
% with an option to export to a sortable html table), or
% tables based on the three different analysis methods:
% indifference-zone ranking, permutation tests, or one-way anova.
%
% Third party code:
% (available from http://www.mathworks.com/matlabcentral/):
% nifti toolbox
% plot3k2.m based on plot3k.m:
%   Ken Garrard, North Carolina State University, 2005
%   Based on plot3c by Uli Theune, University of Alberta
% xticklabel_rotate.m:
%   Author: Denis Gilbert, Ph.D., physical oceanography
%   Maurice Lamontagne Institute, Dept. of Fisheries and Oceans Canada
%   email: gilbertd@dfo-mpo.gc.ca  Web: http://www.qc.dfo-mpo.gc.ca/iml/
%   February 1998; Last revision: 24-Mar-2003
%
% (c) 2008, @rno klein
%--------------------------------------------------------------

%------------------------------------------------------------------------
% Setup: set parameters for input data, evaluation measures, and plotting
%------------------------------------------------------------------------
more off; format compact; 

datapath = '/piece-of-mind/registration_study/';
path(path,'./Utilities/');
path(path,'./Utilities/nifti/');

test_methods = {'FLIRT','AIR','ANIMAL','ART','Demons','FNIRT','ITK','Fluid','Romeo','SICLE','SyN','SPM5_Normalize_direct_8mm','SPM5_Normalize','SPM5_UnifiedSegment','SPM5_DARTEL_pairs'};%,'SPM5_DARTEL'};
test_labels  = {'FLIRT','AIR','ANIMAL','ART','D.Demons','FNIRT','IRTK','JRD-fluid','ROMEO','SICLE','SyN','SPM_N*','SPM_N','SPM_US','SPM_D'};%,'SPM_D'};
%test_labels = test_methods;

category     = 0; % 0 = volume overlap
                  % 1 = surface overlap
                  % 2 = volume similarity
                  % 3 = distance
  overlap_msr  = 1; % for volume overlap: plot measure 1 (TO), 2 (UO), or 3 (MO),
                    % 4 (FP), 5 (FN)
                    
plot_each      = 0; % plot data for each and every registration
plot_boxesROI  = 0; % plot data avg'd across ROIs, then across pairs
                    % -- must have byRegion = 1,
plot_combo     = 0; % plot whole-brain data combining registrations for each method
                    % ex: category=3, byRegion=0
plotMatrix   = 0;     
  plotALLmatrices = 0;
  plotAVGmatrix   = 0;
  plotRanks       = 0; 
    borderfraction= 1/6;
    printRanks    = 0;
byRegion       = 0; % plot data by region: with plot_combo==1 or plotMatrix==1
    orderbysize  = 0;
    plotBrain    = 0;
if plotMatrix==1 || plot_boxesROI==1, byRegion=1; %orderbysize=1; 
end
if plotMatrix==1 && category==2, disp('Need to fix code to plot VS matrices');
end
plot_self      = 0;   % plot resulots of brains registered to themselves
errbar         = 0;   % plot error bars (+-SD) instead of boxplots
lowthresh      = 0.25; % print to stdout when result is below this threshold
lowthresh2     = 0.8; % print to stdout when result is below this threshold
build_table    = 4; % 0 (none), 1 (full table), 
                    % 2 (linear mixed effects model), or 3 (same as 2 for all ROIs)
                    % 4 (ROIxMethod averages) -- must have plotMatrix = 1,
                    %                            plotAVGmatrix = 1
  export_html  = 1;
  if build_table==4, plotMatrix=1; plotAVGmatrix=1; byRegion=1;
  end
  run_anova    = 0;
 
iset_start     = 1;   % start with files beginning with this filestem_tag
iset_end       = 4;   % stop  with files beginning with this filestem_tag

filestem_tags  = {'g';'m';'l';'c'};
tags           = {'MGH10';'CUMC12';'LPBA40';'IBSR18'};
src_filestems1 = {'g1';'g2';'g3';'g4';'g5';'g6';'g7';'g8';'g9';'g10'};
src_filestems2 = {'m1';'m2';'m3';'m4';'m5';'m6';'m7';'m8';'m9';'m10';'m11';'m12'};
src_filestems3 = {'l1';'l2';'l3';'l4';'l5';'l6';'l7';'l8';'l9';'l10';
                  'l11';'l12';'l13';'l14';'l15';'l16';'l17';'l18';'l19';'l20';
                  'l21';'l22';'l23';'l24';'l25';'l26';'l27';'l28';'l29';'l30';
                  'l31';'l32';'l33';'l34';'l35';'l36';'l37';'l38';'l39';'l40'};
src_filestems4 = {'c1';'c2';'c3';'c4';'c5';'c6';'c7';'c8';'c9';'c10';'c11';'c12';'c13';'c14';'c15';'c16';'c17';'c18'};
tgt_filestems1 = src_filestems1;
tgt_filestems2 = src_filestems2;
tgt_filestems3 = src_filestems3;
tgt_filestems4 = src_filestems4;

labels_path    = [datapath 'Labels/'];
in_path        = [datapath 'Evaluations/'];    
out_path       = [datapath 'Evaluations/'];
atlas_path     = [datapath 'Atlases/Atlases_MNIspace/'];
border_path    = [datapath 'Atlases/Atlases_MNIspace_Borders/']; 

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

verbose        = 1;
time0          = clock; 
nmethods       = length(test_methods);
scrsz          = get(0,'ScreenSize');

% Reset source and target files if running ANOVA tests
% on independent subsets of the data
if run_anova==1 
  randtable=0;
  if randtable==1
   PERMpairs = randperm(size(src_filestems1,1));
   PERMpairs = reshape(PERMpairs,size(src_filestems1,1)/2,2);
   src_filestems1 = src_filestems1(PERMpairs(:,1));
   tgt_filestems1 = tgt_filestems1(PERMpairs(:,2));
   PERMpairs = randperm(size(src_filestems2,1));
   PERMpairs = reshape(PERMpairs,size(src_filestems2,1)/2,2);
   src_filestems2 = src_filestems2(PERMpairs(:,1));
   tgt_filestems2 = tgt_filestems2(PERMpairs(:,2));
   PERMpairs = randperm(size(src_filestems3,1));
   PERMpairs = reshape(PERMpairs,size(src_filestems3,1)/2,2);
   src_filestems3 = src_filestems3(PERMpairs(:,1));
   tgt_filestems3 = tgt_filestems3(PERMpairs(:,2));
   PERMpairs = randperm(size(src_filestems4,1));
   PERMpairs = reshape(PERMpairs,size(src_filestems4,1)/2,2);
   src_filestems4 = src_filestems4(PERMpairs(:,1));
   tgt_filestems4 = tgt_filestems4(PERMpairs(:,2));
  else
   src_filestems1 = src_filestems1(1:2:end);
   tgt_filestems1 = tgt_filestems1(2:2:end);
   src_filestems2 = src_filestems2(1:2:end);
   tgt_filestems2 = tgt_filestems2(2:2:end);
   src_filestems3 = src_filestems3(1:2:end);
   tgt_filestems3 = tgt_filestems3(2:2:end);
   src_filestems4 = src_filestems4(1:2:end);
   tgt_filestems4 = tgt_filestems4(2:2:end);
  end
end

% Set number of measures, titles for plots
if category==3, nmeasures = 2;
else            nmeasures = 7;
end

if plot_self==1
     strself = ' (registered to itself)';
else strself = '';
end
                     str0 = '';       strtitle = 'Volume Overlap'; 
if     category==0,  strc = '';       strtitle = 'Volume Overlap';    ylimits = [0.5 1];
elseif category==1,  strc = '_surf';  strtitle = 'Surface Overlap';   ylimits = [0.5 1];
elseif category==2,  strc = '';       strtitle = 'Volume Similarity'; ylimits = [0.5 1];             
elseif category==3,  strc = '_dist';  strtitle = 'Distance Error';    ylimits = [0 10]; 
end

%--------------------------
% Loop through data sources
%--------------------------
for iset = iset_start:iset_end

  s=sprintf('src_filestems = src_filestems%d;', iset); eval(s);
  s=sprintf('tgt_filestems = tgt_filestems%d;', iset); eval(s); 
  npairs       = length(src_filestems) * length(tgt_filestems);
  filestem_tag = char(filestem_tags(iset));
  full_tag     = char(tags(iset));
  
  if build_table==1
       nruns   = size(src_filestems,1);
  else nruns   = (size(src_filestems,1))^2 - size(src_filestems,1);
  end  

%---------------------------
% Prep region-based analysis
%---------------------------
  if byRegion==1 || plotMatrix==1

    % Load label list and filter with labels from first dataset
      isource      = 1;
      itarget      = 1;
      imethod      = 1;
      s=['load ' in_path char(test_methods(imethod)) '/' char(src_filestems(isource)) '_to_' char(tgt_filestems(itarget)) str0 '_eval.mat;'];
      eval(s);
      %   ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN              1x5                40  double              
      %   Labels                                             152x1              1216  double              
      %   label_ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN      152x6              7296  double              
      %   label_nTEST_nTRUE_nAND_nOR                         152x5              6080  double              
      %   nTEST_nTRUE_nAND_nOR                                 1x4                32  double              
      %   ndiff1                                               1x1                 8  double              
      s=sprintf('load %s%s_Labels;',labels_path,char(filestem_tags(iset))); eval(s);
      %   Labels      152x1              1216  double              
      %   Names       152x1             16248  cell                
      I_nTRUE_0                    = find(label_nTEST_nTRUE_nAND_nOR(:,3)>0);
      dummy                        = label_nTEST_nTRUE_nAND_nOR(I_nTRUE_0,1);
      [cLabels,Iintersect,ILabels] = intersect(dummy,Labels);
      cNames                       = Names(ILabels);
      nlabels                      = size(cLabels,1);
    % Use these indices for evaluation
      Ic                           = I_nTRUE_0(Iintersect);
    % Order ROIs by size (NOTE: calculated from first target atlas)
      if orderbysize==1
         [Y,IsortL] = sort(label_nTEST_nTRUE_nAND_nOR(Ic,3),1,'descend');
         cLabels    = cLabels(IsortL);
         cNames     = cNames(IsortL);
         Ic         = Ic(IsortL);
      end     
    % Initialize for 3-D anatomical plots
      if plotBrain==1 || (plotMatrix==1 && plotRanks==1)     
         target_name = char(tgt_filestems(1));
         s=sprintf('TRUE = load_nii(''%s'');',[atlas_path target_name]); eval(s);
         TRUE = double(TRUE.img);
         B = zeros(size(TRUE));
      end      
    % Initialize for matrix plots
      if plotMatrix==1
         if plotRanks==1
           Ranks = zeros(nmethods,nmethods,npairs-length(tgt_filestems),nlabels); %,nmeasures+1);
         end
      end      
        ROIallMethods = zeros(nmethods,npairs-length(tgt_filestems),nlabels,nmeasures+1);
      VSROIallMethods = zeros(nmethods,npairs-length(tgt_filestems),nlabels);
  end
  
% Initialize tables
  if build_table==1
     TOT_ANDoverTRUExMethod = zeros(nruns,nmethods+2); 
  end  
  if build_table==2
   % Method, Brain1, Brain2, Pair, Sign, Value
     fid1=fopen('table_noROI.csv','wt');
     fprintf(fid1, 'Method, Brain1, Brain2, Pair, Sign, Value\n');
  end
  if build_table==3
   % Method, Brain1, Brain2, Pair, Sign, Region, Size, Value
     fid2=fopen('table_ROI.csv','wt');
     fprintf(fid2, 'Method, Brain1, Brain2, Pair, Sign, Region, Size, Value\n');
  end
  tabapp='';
  if build_table==4
   % ROI x Method averages
     if category==0 || category==1
       if     overlap_msr==1, tabapp='TO';
       elseif overlap_msr==2, tabapp='UO';
       elseif overlap_msr==3, tabapp='MO';
       end
       if category==1
          tabapp=[tabapp 'surf'];
       end
     elseif category==2, tabapp='VS';
     elseif category==3, tabapp='DE';
     end
     tabapp=[tabapp '_' full_tag];
     if export_html==1
        fid4=fopen(['table_ROIxMethod_' tabapp '.html'],'wt');
        fprintf(fid4, '<html><head><script type="text/javascript" src="http://www.mindboggle.info/shared/sorttable.js"></script></head><body><span style=font-size:8pt><div style="text-align:left"><table class="sortable" id="sortable_example" style="margin-left:auto">\n');
     else
        fid4=fopen(['table_ROIxMethod_' tabapp '.csv'],'wt');
     end
     stl='';
     for itm=1:nmethods
        if export_html==1
            stl = [stl '<th>' char(test_labels(itm)) '</th>'];
        else
            if itm==nmethods
                 stl = [stl char(test_labels(itm)) ' mean, SD\n'];
            else stl = [stl char(test_labels(itm)) ' mean, SD, '];
            end
        end
     end
     if export_html==1
        fprintf(fid4, ['<tr><th>Region</th>' stl '</tr>']);
     else
        fprintf(fid4, ['Algorithms,' stl]);
     end
  end
 
% Initialize method-total matrix: [method x pairs x measures]
  MTOT = zeros(nmethods,npairs-length(tgt_filestems),nmeasures);

%---------------------
% Loop through methods
%---------------------
  iter_table0 = 1;
  for imethod = 1:nmethods
 
    test_method    = char(test_methods(imethod));
    test_label     = char(test_labels(imethod));
    in_path_method = [in_path test_method '/'];

  % Initialize total and region-based matrices 
    TOT      = zeros(npairs-length(tgt_filestems),nmeasures);
    TOT_self = zeros(       length(tgt_filestems),nmeasures);
    if byRegion==1
       ROI      = zeros(npairs-length(tgt_filestems),nlabels,nmeasures+1);
       ROI_self = zeros(       length(tgt_filestems),nlabels,nmeasures+1);
    end

  %---------------------
  % Loop through targets
  %---------------------
    iter        = 1;
    iter_self   = 1;
    iterNOdep   = 0;
    iterSdep    = 0;
    iterTdep    = 0;
    iterSTdep   = 0;
    ipair_table0 = 1;
    for itarget = 1 : length( tgt_filestems )

      target_name = char(tgt_filestems(itarget));

    %---------------------------------------------------------------
    % For ranking: set the fraction of border voxels for each region
    %---------------------------------------------------------------
      if plotMatrix==1 && plotRanks==1
       %----------------------------
       % Use (LPBA40) target borders
       %----------------------------
         if iset==3
            s=sprintf('load(''%s'');',[border_path target_name '_borders.mat']); eval(s);
            %  Borders      168177x4             5381664  double
            roiborderfractions = zeros(nlabels,2);
            for iLbl=1:nlabels
                IcL     = find(TRUE(:)==cLabels(iLbl));
                sizeIcL = size(IcL,1);
                I       = find(Borders(:,4)==cLabels(iLbl)); 
                if size(I,1)*sizeIcL > 0
                   roiborderfractions(iLbl) = size(I,1)/sizeIcL;
                end
            end
       %--------------------------------------------------------------------------
       % else compute the number of border voxels if the target region were a cube
       %--------------------------------------------------------------------------
         else
            roiborderfractions = zeros(nlabels,2);
            for iLbl=1:nlabels
                IcL            = find(TRUE(:)==cLabels(iLbl));
                sizeIcL        = size(IcL,1);
                cube_edge_len  = nthroot(sizeIcL,3);
                cube_surf      = 6*cube_edge_len^2 - 12*cube_edge_len;
                if cube_surf*sizeIcL > 0
                   roiborderfractions(iLbl) = cube_surf/sizeIcL;
                end
            end
         end
      end

    %---------------------
    % Loop through sources
    %---------------------
      for isource = 1 : length( src_filestems )

        source_name = char(src_filestems(isource));     
         
        isE = exist([in_path_method source_name '_to_' target_name strc '_eval.mat']); 
        if isE>0
          %------------------------------------------------
          % Load source-target registration evaluation file 
          % (generated by evaluate_warps_main.m)
          %------------------------------------------------
            s=['load ' in_path_method source_name '_to_' target_name strc '_eval.mat;'];
            eval(s);
            %   ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN              1x5                40  double              
            %   Labels                                             152x1              1216  double              
            %   label_ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN      152x6              7296  double              
            %   label_nTEST_nTRUE_nAND_nOR                         152x5              6080  double              
            %   nTEST_nTRUE_nAND_nOR                                 1x4                32  double              
            %   ndiff1                                               1x1                 8  double              

          % Change NaNs to zeros -- all but LPBA40 had NaNs in
          % label_ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN
          % and LPBA40 had NaNs in label_meanDist_stdDist
          %   s3 = sum(sum(isnan(label_ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN(:))));
          %   if s3>0, disp([num2str(s3) ' NaNs in s3!']); end
            if category==3
               meanDist_stdDist(isnan(meanDist_stdDist)==1)=0;
               if byRegion==1
                  label_meanDist_stdDist(isnan(label_meanDist_stdDist)==1)=0;
                % Choose non-zero (sorted) labels above   
                  label_meanDist_stdDist = label_meanDist_stdDist(Ic,:);
               end
            else
               nTEST_nTRUE_nAND_nOR(isnan(nTEST_nTRUE_nAND_nOR)==1)=0;
               ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN(isnan(ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN)==1)=0;
            end
            if byRegion==1
              if size(label_nTEST_nTRUE_nAND_nOR,1)>=max(Ic)
               label_nTEST_nTRUE_nAND_nOR(isnan(label_nTEST_nTRUE_nAND_nOR)==1)=0;
               label_ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN(isnan(label_ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN)==1)=0;
             % Choose non-zero (sorted) labels above 
               label_nTEST_nTRUE_nAND_nOR = label_nTEST_nTRUE_nAND_nOR(Ic,:);
               label_ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN = label_ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN(Ic,:);
               lentest=1;
              else disp('label_nTEST_nTRUE_nAND_nOR and max(Ic) mismatch:');
                   disp([isource itarget]);
                   lentest=0;
              end
            end
            
          %----------------------
          % Case: target = source
          %----------------------
            if itarget==isource
                
             %--------------------------------------
             % Build total and region-based matrices
             %--------------------------------------
               if byRegion==1
                  if category==3
                       ROI_self(iter_self,:,:) = label_meanDist_stdDist;
                  %elseif category==2, disp('FIX');
                  else
                      if lentest==1
                        ROI_self(iter_self,:,:) = [label_ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN ...
                                                  label_nTEST_nTRUE_nAND_nOR(:,2:3)];
                      end
                  end
               end
               if category==3
                    TOT_self(iter_self,:) = meanDist_stdDist;
               else TOT_self(iter_self,:) = [ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN ...
                                             nTEST_nTRUE_nAND_nOR(1:2)];
               end
               iter_self = iter_self + 1;
               
          %-----------------------
          % Case: target ~= source
          %-----------------------               
            else
                
             %-------------
             % Build tables
             %-------------
               if build_table==1
                  TOT_ANDoverTRUExMethod(iter,1:2) = [isource itarget];
                  if     category==0 || category==1
                     %TOT_ANDoverTRUExMethod(iter,imethod+2) = ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN(overlap_msr);
                     TOT_ANDoverTRUExMethod(iter,imethod+2) = mean(label_ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN(:,overlap_msr+1));
                  elseif category==2, disp('FIX');
                  elseif category==3 
                     TOT_ANDoverTRUExMethod(iter,imethod+2) = meanDist_stdDist(1);
                  end                      
               end              
			   if build_table==2
                  if     itarget>isource, signofpair = '1';
                  elseif itarget<isource, signofpair = '-1';
                  end
                % Method, Brain1, Brain2, Pair, Sign, Value
                  fprintf(fid1, '"%d", "%d", "%d", "%d,%d", "%s", "%f", "%f", "%f", "%f", "%f"\n',...
                      imethod, isource, itarget, isource, itarget, signofpair,...
                      ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN); 
               end  
               
             %-------------------
             % Build total matrix
             %-------------------
               if category==3
                    TOT(iter,:) = meanDist_stdDist;
               elseif category==2
                    TOT(iter,:) = [ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN ...
                                   nTEST_nTRUE_nAND_nOR(1:2)];
                    vs    = 1 - abs((TOT(iter,6)-TOT(iter,7))) ./ (TOT(iter,6)+TOT(iter,7));
                    if vs<lowthresh2
                       if verbose > 0
                          s=sprintf('disp(''%s:  LOW VS for %s_to_%s: %1.2f'');',test_method,source_name,target_name,vs); eval(s);
                       end
                    end                   
               else TOT(iter,:) = [ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN ...
                                   nTEST_nTRUE_nAND_nOR(1:2)];
                    if ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN(1)<lowthresh
                       if verbose > 0
                          s=sprintf('disp(''%s:  LOW TO for %s_to_%s: %1.2f'');',test_method,source_name,target_name,ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN(1)); eval(s);
                       end
                    end
               end
               
             %--------------------------
             % Build region-based matrix
             %--------------------------              
               if byRegion==1
                  if category==3
                       ROI(iter,:,:) = label_meanDist_stdDist;
                  else
                      if lentest==1
                        ROI(iter,:,:) = [label_ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN ...
                                        label_nTEST_nTRUE_nAND_nOR(:,2:3)];                             
                      end
                      %label_ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN
                 %-------------------------
                 % Build region-based table
                 %-------------------------              
                   if build_table==3
                     if     itarget>isource, signofpair = '1';
                     elseif itarget<isource, signofpair = '-1';
                     end
                     fprintf(fid2, '"%d", "%d", "%d", "%d,%d", "%s", "%d", "%d", "%f"\n',...
                             imethod, isource, itarget, isource,itarget, signofpair, ...
                             0,...
                             nTEST_nTRUE_nAND_nOR(1),... 
                             ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN(1)); 
                   % Method, Brain1, Brain2, Pair, Sign, Region, Size, Value
                     for iROI=1:nlabels
                       fprintf(fid2, '"%d", "%d", "%d", "%d,%d", "%s", "%d", "%d", "%f"\n',...
                              imethod, isource, itarget, isource,itarget, signofpair, ...
                              label_nTEST_nTRUE_nAND_nOR(iROI,1),...
                              label_nTEST_nTRUE_nAND_nOR(iROI,3),... 
                              label_ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN(iROI,2)); 
                     end
                   end
                  end
               end
               iter = iter + 1;
            end
        end
      end   % for isource   
    end     % for itarget

  %-----------
  % Save table
  %-----------
	if build_table==1
       if run_anova==0
  	      %s=sprintf('save table%d.txt TOT_ANDoverTRUExMethod -ascii -double -tabs',nruns); eval(s);
  	      s=sprintf('save table%d_avgTO.txt TOT_ANDoverTRUExMethod -ascii -double -tabs',nruns); eval(s);
       end
    end
    
  %-------------
  % Set up plots
  %-------------
    if byRegion==1
       ROIallMethods(imethod,:,:,:) = ROI;
    end

    if plot_each==1 || (plot_combo==1 && byRegion==1)
         titlec = [strtitle strself ' for ' test_label ' (' full_tag ')'];
    else titlec = [strtitle strself ' (' full_tag ')'];
    end
    if category==0    
         ylabelc       = 'Overlap'; %Volume (1mm^3 voxels in MNI152 space)';  
    elseif category==1
         ylabelc       = 'Overlap'; %Area (1mm^3 voxels in MNI152 space)';
    elseif category==2            
         ylabelc       = 'Volume Similarity'; %Volume (1mm^3 voxels in MNI152 space)';
         if plot_each==1
              titlec = [strtitle strself ' for ' test_label ' (' full_tag ')'];
         else titlec = [strtitle strself ' (' full_tag ')'];
         end
    elseif category==3
         ylabelc       = 'Distance (mm in MNI152 space)';
    end
    
    if plot_self==1,  PTOT = TOT_self;  if byRegion==1, PROI = ROI_self; end
    else              PTOT = TOT;       if byRegion==1, PROI = ROI;      end
    end

  %----------------------------
  % Calculate volume similarity
  %----------------------------
    if category==2,    VS    = 1 - abs(PTOT(:,6)-PTOT(:,7)) ./ (PTOT(:,6)+PTOT(:,7));
       if byRegion==1, VSroi = 1 - abs(PROI(:,:,7)-PROI(:,:,8)) ./ (PROI(:,:,7)+PROI(:,:,8));
                       VSROIallMethods(imethod,:,:) = VSroi;
       end
       if mean(VS)<lowthresh2
          if verbose > 0
             s=sprintf('disp(''%s:  LOW VS for %s_to_%s: %1.2f'');',test_method,source_name,target_name,mean(VS)); eval(s);
          end
       end
    end
    
  %-----------
  % PLOT EACH
  %-----------  
    if plot_each > 0 
        
      figure(iset+1000+category*10); 
      subplot(nmethods+1,1,imethod), 
      if category==2
         %bar(PTOT(:,[6 7]),'grouped');
          bar(VS);
          axis([1 size(PTOT,1) ylimits]);
      elseif category==3
          bar(PTOT(:,1));
          axis([1 size(PTOT,1) ylimits]);
      else
          plot(PTOT(:,1:3));
          if imethod==nmethods
             legend('AND/TRUE','AND/OR','AND/MEAN'); 
          end
          axis([1 size(PTOT,1) ylimits]);    
      end
      title(titlec);
      
    end

  %--------------------
  % PLOT COMBINED DATA
  %--------------------   
    if plot_combo==1
     %-------------
     % Total graphs
     %-------------
       if byRegion==0

          figure(iset+250+category*10+byRegion); hold on

          if errbar==0
             if category==2, MTOT(imethod,:,1) = VS;    
             else
                 MTOT(imethod,:,:) = PTOT;
                 if imethod==1
                   disp(full_tag);
                 end
                 a=[char(test_labels(imethod)) ' ' num2str(mean(PTOT(:,overlap_msr))) ' ' num2str(std(PTOT(:,overlap_msr)))];
                 disp(a);
             end
             if imethod==nmethods
                axis([0 nmethods+1 ylimits]);
                set(gca,'XTick',1:nmethods,'YLim',ylimits); %,'XTickLabel',test_labels)
                if nmethods>1
                    xticklabel_rotate(1:nmethods,60,test_labels,'interpreter','none')
                end
             end
          else
             if     category==2
                    errorbar(imethod,mean(VS),std(VS,0,1),'b.');
             elseif category==3
                    Inotnan = find(isnan(PROI(:,iPROI,2))==0);
                    errorbar(imethod,mean(PTOT(Inotnan,1),1),std(PTOT(Inotnan,1),0,1),'b.');
             else   errorbar(imethod,mean(PTOT(:,overlap_msr)),std(PTOT(:,overlap_msr)),'b.');
             end
             title(titlec);
             if imethod==nmethods
                ylabel(ylabelc);
                axis([0 nmethods+1 ylimits]);
                set(gca,'XTick',1:nmethods,'XTickLabel',test_labels,'YLim',ylimits)
             end
          end
          
     %--------------------
     % Region-based graphs
     %--------------------
       else

          figure(iset+category*10+byRegion);  
          subplot(nmethods,1,imethod); hold on

          for iPROI=1:size(PROI,2)
         
              if category==2
                  errorbar(iPROI,mean(VSroi(:,iPROI,1)),std(VSroi(:,iPROI,1)),'b.');                   
                  if plotBrain==1
                     I=find(TRUE==PROI(1,iPROI,1));
                     if size(I,1)>0
                       B(I) = mean(VSroi(:,iPROI,1));
                     end                    
                  end
              elseif category==3
                  Inotnan = find( isnan(PROI(:,iPROI,2)) + isinf(PROI(:,iPROI,2)) == 0 );
                  errorbar(iPROI,mean(PROI(Inotnan,iPROI,2)),std(PROI(Inotnan,iPROI,2)),'b.');
                  if plotBrain==1
                     I=find(TRUE==PROI(1,iPROI,1));
                     if size(I,1)>0
                        B(I) = mean(PROI(Inotnan,iPROI,2));
                       %B(I) = (mean(PROI(Inotnan,iPROI,2)))/max(PROI(Inotnan,iPROI,2));
                     end                    
                  end
              else
                  errorbar(iPROI,mean(PROI(:,iPROI,overlap_msr+1)),std(PROI(:,iPROI,overlap_msr+1)),'b.');
                  set(gca,'YLim',ylimits);
                  if plotBrain==1
                     I=find(TRUE==PROI(1,iPROI,1));
                     if size(I,1)>0
                       B(I) = mean(PROI(:,iPROI,2));
                     end                    
                  end 
              end             
              title(titlec);
              axis([0 nlabels+1 ylimits]);
              if imethod==nmethods
                 if iPROI==nlabels
                    set(gca,'XTick',1:nlabels,'XTickLabel',cNames)
                    xticklabel_rotate(1:nlabels,90,cNames,'interpreter','none')
                 end
              end   
          end
          
        %-----------------
        % Anatomical plots
        %-----------------
          if plotBrain==1
              mapjet = colormap('jet');            
              if category~=3
                 mapjet = mapjet(:,end:-1:1);
              end

              figure(iset+category*10+byRegion+imethod+plotBrain*1000); hold on
              unique(B(:))
              I = find(B~=0); % & B<1);
              [x,y,z] = ind2sub(size(TRUE),I);
              plot3k2([x,y,z],B(I),ylimits(1),ylimits(2),mapjet);
              axis equal; axis off
              title(titlec);
              view(-100,30)
              print_tiffs=0;
              if print_tiffs==1
                  if category==3
                    s=sprintf('print -dtiff -r600 Brain_DE_%s_%s.tiff',...
                    full_tag,char(test_labels(imethod))); eval(s); 
                  elseif category==2
                    s=sprintf('print -dtiff -r600 Brain_VS_%s_%s.tiff',...
                    full_tag,char(test_labels(imethod))); eval(s); 
                  elseif category==1
                    if overlap_msr==1
                      s=sprintf('print -dtiff -r600 Brain_TOsurf_%s_%s.tiff',...
                      full_tag,char(test_labels(imethod))); eval(s); 
                    elseif overlap_msr==2
                      s=sprintf('print -dtiff -r600 Brain_U0surf_%s_%s.tiff',...
                      full_tag,char(test_labels(imethod))); eval(s); 
                    elseif overlap_msr==3
                      s=sprintf('print -dtiff -r600 Brain_MOsurf_%s_%s.tiff',...
                      full_tag,char(test_labels(imethod))); eval(s); 
                    end
                  else
                    if overlap_msr==1
                      s=sprintf('print -dtiff -r600 Brain_TO_%s_%s.tiff',...
                      full_tag,char(test_labels(imethod))); eval(s); 
                    elseif overlap_msr==2
                      s=sprintf('print -dtiff -r600 Brain_U0_%s_%s.tiff',...
                      full_tag,char(test_labels(imethod))); eval(s); 
                    elseif overlap_msr==3
                      s=sprintf('print -dtiff -r600 Brain_MO_%s_%s.tiff',...
                      full_tag,char(test_labels(imethod))); eval(s); 
                    end
                  end
              end
          end
       end
    end

  %----------------------------------------
  % Matrix (color table) plots: everything!
  %----------------------------------------
    if plotMatrix==1 && plotALLmatrices==1
         figure(10+iset)
         subplot(1,nmethods,imethod)
         imagesc(squeeze(ROIallMethods(imethod,:,:,overlap_msr+1))',[0 1]); 
         mapjet = colormap('jet');            
         if category~=3
            mapjet = mapjet(:,end:-1:1);
         end           
         title([test_label ': ROIs x Runs (' full_tag ')']);
         colormap(mapjet); %colorbar
         axis off
    end

  end  % for imethod = 1:nmethods 

%----------------
% Total box plots: averages
%---------------- 
  if byRegion==1
    if plot_boxesROI==1
          ylimits_box = [0 1];
          if category==0
            meanboxes = squeeze(mean(ROIallMethods(:,:,:,overlap_msr+1),3));
          elseif category==3
            meanboxes = squeeze(mean(ROIallMethods(:,:,:,2),3));
            ylimits_box = [0 10]; axis ij
          elseif category==2
            meanboxes = squeeze(mean(VSROIallMethods(:,:,:),3));
          end
          figure(10000+100*category+10*iset+1)
          boxplotC(meanboxes',0,'+',1,1.5,'k',0,1.25);
          title(titlec);
          xlabel('');
          ylabel(ylabelc);
          if imethod==nmethods
            axis([0 nmethods+1 ylimits_box]);
            set(gca,'XTick',1:nmethods,'YLim',ylimits_box); %,'XTickLabel',test_labels)
            if nmethods>1
                xticklabel_rotate(1:nmethods,60,test_labels,'interpreter','none')
            end
          end
    end
  end

%----------------
% Total box plots: whole-brain
%----------------  
  if plot_combo==1 && byRegion==0 && errbar==0
      if category==2
         boxplot(MTOT(:,:,1)','notch','on');%,'labels',test_labels);
      elseif category==3
          boxplotC(MTOT(:,:,1)',1,[],1,'k','k',false,2);
          %boxplot([MTOT(:,:,1)'],'notch','on');%,'labels',test_labels);
      else
         fh=boxplot(MTOT(:,:,overlap_msr)','notch','on','symbol','k+','colors','k'); hold on %,'labels',test_labels); hold on
      end             
      title(titlec);
      xlabel('');
      ylabel(ylabelc);
      axis([0 nmethods+1 ylimits]);
      %set(gca,'XTick',[1:nmethods],'XTickLabel',test_labels)
  end
    
%---------------------------
% Matrix (color table) plots
%---------------------------
  if plotMatrix==1

  %------------------
  % Rank matrix plots
  %------------------
    if plotRanks==1 

        for ipairr = 1:size(ROI,1)
            for iroir = 1:nlabels
              if     category==0, ievalr = overlap_msr;
              elseif category==1, ievalr = overlap_msr; 
              elseif category==2, error('FIX'); 
              elseif category==3, ievalr = 1;
              end                            
              R     = ROIallMethods(:,ipairr,iroir,ievalr+1);
              delta = borderfraction * roiborderfractions(iroir);
              R     = rank_indifferencezone(R,delta);
              Ranks(:,:,ipairr,iroir) = R;              
            end
        end

        saveRanks=0;
        if saveRanks==1
          if category==3
            s=sprintf('save RankMatrix_ROI_DE_%s.mat',full_tag); eval(s);
          elseif category==2
            s=sprintf('save RankMatrix_ROI_VS_%s.mat',full_tag); eval(s);
          elseif category==1
            if overlap_msr==1
              s=sprintf('save RankMatrix_ROI_TOsurf_%s.mat',full_tag); eval(s);
            elseif overlap_msr==2
              s=sprintf('save RankMatrix_ROI_UOsurf_%s.mat',full_tag); eval(s);
            elseif overlap_msr==3
              s=sprintf('save RankMatrix_ROI_MOsurf_%s.mat',full_tag); eval(s);
            end
          else
            if overlap_msr==1
              s=sprintf('save RankMatrix_ROI_TO_%s.mat',full_tag); eval(s);
            elseif overlap_msr==2
              s=sprintf('save RankMatrix_ROI_UO_%s.mat',full_tag); eval(s);
            elseif overlap_msr==3
              s=sprintf('save RankMatrix_ROI_MO_%s.mat',full_tag); eval(s);
            end
          end
        end

        mean_acrosspairs = squeeze(mean(squeeze(mean(Ranks,3)),2));
        mean_acrossROIs  = squeeze(mean(squeeze(mean(Ranks,4)),2));

        disp('avg across ROIs (after avg across pairs):');
        meanRanks_acrossROIs = mean(mean_acrosspairs,2);
        stdRanks_acrossROIs  = std(mean_acrosspairs,[],2);
        [S,IsortR] = sort(meanRanks_acrossROIs,'descend');
        disp([char(test_labels(IsortR)) num2str(meanRanks_acrossROIs(IsortR)) num2str(stdRanks_acrossROIs(IsortR))]);

        disp('avg across pairs (after avg across ROIs):');
        meanRanks_acrosspairs = mean(mean_acrossROIs,2);
        stdRanks_acrosspairs  = std(mean_acrossROIs,[],2);
        [S,IsortR] = sort(meanRanks_acrosspairs,'descend');
        disp([char(test_labels(IsortR)) num2str(meanRanks_acrosspairs(IsortR)) num2str(stdRanks_acrosspairs(IsortR))]);

        figure('Position',[1 scrsz(4) scrsz(3)*0.2 scrsz(4)]);
        imagesc(mean_acrosspairs');
        mapjet = colormap('jet');            
        if category~=3
           mapjet = mapjet(:,end:-1:1);
        end
        colormap(mapjet)
        h=title({'Ranked methods: overlap';[full_tag ' (' num2str(size(Ranks,3)) 'pairs)']});
        set(h,'FontSize',12,'FontWeight','Bold');
        set(gca,'YTick',1:nlabels,'YTickLabel',cNames,'FontSize',5,'FontWeight','Bold')
        set(gca,'XTick',1:nmethods,'XTickLabel',test_labels,'FontSize',10,'FontWeight','Bold')
        set(gca,'TickLength',[0 0]);
        xticklabel_rotate(1:nmethods,60,test_labels,'interpreter','none')
        colorbar       

        if printRanks==1
          if category==3
            s=sprintf('print -dtiff -r600 RankMatrix_ROI_DE_%s.tiff',full_tag); eval(s); 
            s=sprintf('saveas(gcf,''RankMatrix_ROI_DE_%s.fig'');',full_tag); eval(s); 
          elseif category==2
            s=sprintf('print -dtiff -r600 RankMatrix_ROI_VS_%s.tiff',full_tag); eval(s); 
            s=sprintf('saveas(gcf,''RankMatrix_ROI_VS_%s.fig'');',full_tag); eval(s); 
          elseif category==1
            if overlap_msr==1
              s=sprintf('print -dtiff -r600 RankMatrix_ROI_TOsurf_%s.tiff',full_tag); eval(s); 
              s=sprintf('saveas(gcf,''RankMatrix_ROI_TOsurf_%s.fig'');',full_tag); eval(s); 
            elseif overlap_msr==2
              s=sprintf('print -dtiff -r600 RankMatrix_ROI_UOsurf_%s.tiff',full_tag); eval(s); 
              s=sprintf('saveas(gcf,''RankMatrix_ROI_UOsurf_%s.fig'');',full_tag); eval(s); 
            elseif overlap_msr==3
              s=sprintf('print -dtiff -r600 RankMatrix_ROI_MOsurf_%s.tiff',full_tag); eval(s); 
              s=sprintf('saveas(gcf,''RankMatrix_ROI_MOsurf_%s.fig'');',full_tag); eval(s); 
            end
          else
            if overlap_msr==1
              s=sprintf('print -dtiff -r600 RankMatrix_ROI_TO_%s.tiff',full_tag); eval(s); 
              s=sprintf('saveas(gcf,''RankMatrix_ROI_TO_%s.fig'');',full_tag); eval(s); 
            elseif overlap_msr==2
              s=sprintf('print -dtiff -r600 RankMatrix_ROI_UO_%s.tiff',full_tag); eval(s); 
              s=sprintf('saveas(gcf,''RankMatrix_ROI_UO_%s.fig'');',full_tag); eval(s); 
            elseif overlap_msr==3
              s=sprintf('print -dtiff -r600 RankMatrix_ROI_MO_%s.tiff',full_tag); eval(s); 
              s=sprintf('saveas(gcf,''RankMatrix_ROI_MO_%s.fig'');',full_tag); eval(s); 
            end
          end
        end 
    end
  %----------------------
  % Plot ALL/AVG matrices
  %----------------------
    if plotAVGmatrix==1   
      %-----------------------------------
      % Build ROIxMethod table of averages
      %-----------------------------------    
        if build_table==4
            meanROIs = 100*squeeze(mean(ROIallMethods,2));
            stdROIs  = 100*squeeze(std(ROIallMethods,[],2));
            for itl=1:nlabels
              st2='';
              for itm=1:nmethods
                s=sprintf('mR = %1.2f;',meanROIs(itm,itl,overlap_msr+1)); eval(s);
                s=sprintf('sR = %1.2f;',stdROIs(itm,itl,overlap_msr+1)); eval(s);
                if export_html==1
                    mRnum = num2str(round(mR*255/100));
                    mRstr = [mRnum ',' mRnum ',' mRnum];
                    st2 = [st2 '<td style=background-color:rgb(' mRstr ');color:rgb(0,0,0);padding-left:10px>' num2str(round(mR)) '</td>'];
                else
                    if itm==nmethods
                         st2 = [st2 num2str(mR) ', ' num2str(sR) ' '];
                    else st2 = [st2 num2str(mR) ', ' num2str(sR) ', '];
                    end
                end
              end
              if export_html==1
                 fprintf(fid4, '<tr><td>%s</td>%s</tr>', char(cNames(itl)),st2);
              else
                 fprintf(fid4, '"%s", %s\n', char(cNames(itl)),st2);
              end
            end
            if export_html==1
               fprintf(fid4, '</table></div></span></body></html>');
            end                 
        end
        
      %------------------
      % Plot color tables
      %------------------
        plot_others = 0;
        figure
        if plot_others==1
          subplot(1,5,1)
        end
        meanROIs = squeeze(mean(ROIallMethods,2));
        imagesc(squeeze(meanROIs(:,:,overlap_msr+1))'); 
        mapjet = colormap('jet');            
        if category==3
           mapjet = mapjet(:,end:-1:1);
        end           
        colormap(mapjet);
        title(['AVERAGE AND/TRUE: Regions x Methods (' full_tag ')']);
        set(gca,'YTick',1:nlabels,'YTickLabel',cNames)
        set(gca,'XTick',1:nmethods,'XTickLabel',test_labels)
        set(gca,'TickLength',[0 0]);
        xticklabel_rotate(1:nmethods,90,test_labels,'interpreter','none')
        if plot_others==1
            subplot(1,5,2)
            maxROIs = squeeze(max(ROIallMethods,[],2));
            imagesc(squeeze(maxROIs(:,:,overlap_msr+1))'); 
            mapjet = colormap('jet');            
            if category==3
               mapjet = mapjet(:,end:-1:1);
            end           
            colormap(mapjet);
            title(['MAXIMUM AND/TRUE: Regions x Methods (' full_tag ')']);
            set(gca,'XTick',1:nmethods,'XTickLabel',test_labels)
            set(gca,'TickLength',[0 0]);
            xticklabel_rotate(1:nmethods,90,test_labels,'interpreter','none')
            axis off
            subplot(1,5,3)
            minROIs = squeeze(min(ROIallMethods,[],2));
            imagesc(squeeze(minROIs(:,:,overlap_msr+1))'); 
            mapjet = colormap('jet');            
            if category==3
               mapjet = mapjet(:,end:-1:1);
            end           
            colormap(mapjet);
            title(['MINIMUM AND/TRUE: Regions x Methods (' full_tag ')']);
            set(gca,'XTick',1:nmethods,'XTickLabel',test_labels)
            set(gca,'TickLength',[0 0]);
            xticklabel_rotate(1:nmethods,90,test_labels,'interpreter','none')
            axis off
            subplot(1,5,4)
            stdROIs = squeeze(std(ROIallMethods,[],2));
            imagesc(squeeze(stdROIs(:,:,overlap_msr+1))'); 
            mapjet = colormap('jet');            
            if category==3
               mapjet = mapjet(:,end:-1:1);
            end           
            colormap(mapjet);
            title(['STANDARD DEVIATION AND/TRUE: Regions x Methods (' full_tag ')']);
            set(gca,'XTick',1:nmethods,'XTickLabel',test_labels)
            set(gca,'TickLength',[0 0]);
            xticklabel_rotate(1:nmethods,90,test_labels,'interpreter','none')
            axis off
            subplot(1,5,5)
            andtrueXregsXrois = squeeze(ROIallMethods(:,:,:,2));
            ROIbestMethod     = zeros(size(andtrueXregsXrois));
            for ireg = 1 : size(ROI,1)
                for iroi = 1 : size(ROI,2)          
                    [s,IsortM] = sort(squeeze(andtrueXregsXrois(:,ireg,iroi)));
                    [s,IsortM] = sort(IsortM);
                    ROIbestMethod(:,ireg,iroi) = IsortM;
                end
            end
            imagesc(squeeze(mean(ROIbestMethod,2))'); 
            mapjet = colormap('jet');            
            if category==3
               mapjet = mapjet(:,end:-1:1);
            end           
            colormap(mapjet);
            title(['AVERAGE RANK AND/TRUE: Regions x Methods (' full_tag ')']);
            set(gca,'XTick',1:nmethods,'XTickLabel',test_labels)
            set(gca,'TickLength',[0 0]);
            xticklabel_rotate(1:nmethods,90,test_labels,'interpreter','none')
            axis off
        end
    end 
  end
  
%-----------------------
% ANOVA tables and plots
%-----------------------
  if build_table==1 && run_anova==1
    [p,a,s]=anova1(TOT_ANDoverTRUExMethod,test_labels);
    figure; [c,m,h,nms]=multcompare(s,'estimate','anova1','ctype','bonferroni');
  end

  if build_table==4
     fclose(fid4);
  end
  
end  % for iset = iset_start:iset_end

disp(['Completed in ' num2str(etime(clock,time0)/60) ' minutes']);