function evaluate_distances_main();
% This program was used to generate the label boundary distance errors
% 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 for evaluating automated vs. manual brain label distances
% calls evaluate_distances.m.
%
% (c) 2008, @rno klein
%--------------------------------------------------------------

%---------------------------------------------------
% Setup: load parameters and create output directory
%---------------------------------------------------
 more off; format compact; 
 datapath = '/piece-of-mind/registration_study/';
 
%reg_method  = 'FLIRT'
%reg_method  = 'ART'
%reg_method  = 'ITK'
%reg_method  = 'AIR'
%reg_method  = 'ANIMAL' 
%reg_method  = 'Demons' 
%reg_method  = 'Fluid'
%reg_method  = 'FNIRT' 
reg_method  = 'SyN'
%reg_method  = 'SPM5_Normalize_direct_8mm' % 183x219x183 -> 182x218x182 
%reg_method  = 'SICLE'  
%reg_method  = 'SPM5_Normalize'
%reg_method  = 'Romeo'
%reg_method  = 'Pasha'
%reg_method  = 'HAMMER'
%reg_method  = 'SPM5_UnifiedSegment'
%reg_method  = 'SPM5_DARTEL_pairs' 
 
 override_flip    = 0; % Override flip for 'l' filestem
 flip_true        = 0; 
 if strcmp(reg_method,'SyN')==1 || strcmp(reg_method,'Romeo')==1 || strcmp(reg_method,'SPM5_Normalize_direct_8mm')==1 || strcmp(reg_method,'SPM5_UnifiedSegment')==1 || strcmp(reg_method,'SPM5_DARTEL')==1  || strcmp(reg_method,'SPM5_DARTEL_pairs')==1 
     flip_true    = 1; 
 end
 
 verbose          = 2; % 1=filenames; 2=filenames & values
 save_results     = 0;
 plot_test        = 0;
 
 iset_start       = 3;
 iset_end         = 3;
 
 atlas_path  = [datapath 'Atlases/Atlases_MNIspace_Borders/'];     
 test_path   = [datapath 'Evaluations/'             reg_method '/'];
 out_path    = [datapath 'Evaluations/'             reg_method '/'];
 labels_path = [datapath 'Labels/'];
 
 filestem_tags  = {'g';'m';'l';'c'};
 
 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;
  
 dims1 = [182 218 182];
 dims2 = [182 218 182];
 dims4 = [182 218 182];
 if strcmp(reg_method,'SPM5_UnifiedSegment_head')==1
    dims3 = [256 124 256];
 else
    dims3 = [181 217 181];
 end
 
 time0 = clock; 

 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);
   s=sprintf('dims          = dims%d;',          iset); eval(s);
   s=sprintf('load %s%s_Labels;',labels_path,char(filestem_tags(iset))); eval(s);

   Z = zeros(dims);
   
 %-----------------------------------------------
 % Loop through sources and targets (and atlases)
 %-----------------------------------------------
   for itarget = 1 : length( tgt_filestems )

     target_name = char(tgt_filestems(itarget));

     s=sprintf('load(''%s'');',[atlas_path target_name '_borders.mat']); eval(s);
     %     bfactor = 0.6;
     %     Borders(:,2) = bfactor * Borders(:,2); 
     %     Borders = round(Borders);
     %     [min(Borders,[],1); max(Borders,[],1)]
     I = sub2ind(dims,Borders(:,1),Borders(:,2),Borders(:,3));
     TRUE = Z;
     TRUE(I) = Borders(:,4);

     if flip_true==1
       TRUE = TRUE(end:-1:1,:,:);
     end
          
     for isource = 1 : length( src_filestems )

       source_name =  char(src_filestems(isource));                
       
       %--------------
       % Image volumes
       %--------------
         if verbose > 0
            s=sprintf('disp(''%s_to_%s'');',source_name,target_name); eval(s);
         end

         s=sprintf('load(''%s'');',[test_path source_name '_to_' target_name '_borders.mat']); eval(s);
 
         I = sub2ind(dims,Borders(:,1),Borders(:,2),Borders(:,3));
     
         TEST = Z;
         TEST(I) = Borders(:,4);
        
       %----------
       % Test plot
       %----------
         if plot_test==1
           for ii=1:125
               subplot(1,2,2), imagesc(TEST(:,:,ii)); axis equal; axis off; 
               subplot(1,2,1), imagesc(TRUE(:,:,ii)); axis equal; axis off; 
               colormap(gray); drawnow; pause(0);
           end
         end

       %-----------------------------------------
       % Evaluate labels assigned by the template
       %-----------------------------------------         
         fixit=0;
         if fixit==1
            eval(['load ' out_path source_name '_to_' target_name '_dist_eval.mat']);
            meanDist_stdDist = evaluate_distances_whole( TEST, TRUE, Labels, verbose );
         else
           [label_meanDist_stdDist, meanDist_stdDist] = evaluate_distances( TEST, TRUE, Labels, verbose );
         end
         
         if save_results==1
            eval(['save ' out_path source_name '_to_' target_name '_dist_eval.mat label_meanDist_stdDist meanDist_stdDist;']);
         end
     end
   end
 end

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