function evaluate_warps_main();
% This program was used to generate the label overlap measures
% 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 overlaps
% calls evaluate_warps.m.
%
% Third party code:
% (available from http://www.mathworks.com/matlabcentral/):
% nifti toolbox
%
% (c) 2008, @rno klein
%--------------------------------------------------------------

%---------------------------------------------------
% Setup: load parameters and create output directory
%---------------------------------------------------
 %reg_method  = 'FLIRT'
 %reg_method  = 'ITK'
 %reg_method  = 'Romeo'
 %reg_method  = 'AIR'
 %reg_method  = 'ANIMAL' 
 %reg_method  = 'Demons' 
 %reg_method  = 'Fluid' 
 %reg_method  = 'FNIRT'
 %reg_method  = 'ART'
 %reg_method  = 'SPM5_Normalize_direct_8mm' % 183x219x183 -> 182x218x182
 %reg_method  = 'SICLE'  
 %reg_method  = 'Pasha'
 %reg_method  = 'SPM5_Normalize'
 %reg_method  = 'SPM5_UnifiedSegment'
 %reg_method  = 'SPM5_DARTEL_pairs' 
 reg_method  = 'SyN'
 
 datapath       = '/piece-of-mind/studies/registration_study/';
 test_path      = [datapath 'Transformed_Atlases/' reg_method '/'];
 out_path       = [datapath 'Evaluations/'         reg_method '/'];
 labels_path    = [datapath 'Labels/'];
 % whole-head data: SPM5_*_head
 % atlas_path   = [datapath 'Atlases/Atlases/']; % for SPM5_*_head (native space) 
 % atlas_path   = [datapath 'Atlases/Atlases_MNIspace_RL_DARTELeval/'];
 % atlas_path   = [datapath 'Atlases/Atlases_MNIspace_RL_DARTELpairseval/'];
 atlas_path     = [datapath 'Atlases/Atlases_MNIspace/'];
 % surface_path = [datapath 'Surfaces/Surfaces/'];
 surface_path   = [datapath 'Surfaces/Surfaces_MNIspace/'];
 path(path, '/piece-of-mind/studies/registration_study/Programs/nifti');
 
 save_results     = 0;
 plot_test        = 0;
 
 iset_start       = 4;
 iset_end         = 4;
 flip_true        = 1;  % flip relevant ('l' filestem) data 
 filestem_tags  = {'g';'m';'c';'l'};
 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 = {'c1';'c2';'c3';'c4';'c5';'c6';'c7';'c8';'c9';
                   'c10';'c11';'c12';'c13';'c14';'c15';'c16';'c17';'c18'};
 src_filestems4 = {'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'};
 tgt_filestems1 = src_filestems1;
 tgt_filestems2 = src_filestems2;
 tgt_filestems3 = src_filestems3;
 tgt_filestems4 = src_filestems4; 
 %src_filestems4 = {'l1'};
 %tgt_filestems4 = {'l1';'l11';'l21';'l31'};
 
%----------------------  

 more off; format compact; 
 verbose      = 2; % 1=filenames; 2=filenames & values   
 mask_surface = 0;
   flip_mask  = 1;

 eval(['isE = exist('''  out_path ''',''dir'');']);
 if   isE == 0,  eval(['! mkdir ' out_path ';']);
 end
 if plot_test==1,  figure
 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('load %s%s_Labels;',labels_path,char(filestem_tags(iset))); eval(s);

 %-------------------------------
 % Loop through targets (atlases)
 %-------------------------------
   for itarget = 1 : length( tgt_filestems )

     target_name = char(tgt_filestems(itarget));

     s=sprintf('TRUE = load_nii(''%s'');',[atlas_path target_name]); eval(s);
     TRUE = double(TRUE.img);
     %loaddims  = '[256 124 256]';
     %s=sprintf('fid = fopen(''%s'',''r'',''%s'');',[atlas_path target_name '.img'],'ieee-be'); 
     %eval(s);
     %TRUE = fread(fid,'int16');
     %s=sprintf('TRUE = reshape(TRUE,%s);',loaddims); eval(s);

     if flip_true==1
       TRUE = TRUE(end:-1:1,:,:);
     end
     
     if mask_surface == 1
       disp('Masking surface...');
       s=sprintf('surfMASK = load_nii(''%s'');',[surface_path target_name]); eval(s);
       surfMASK = surfMASK.img>0;
       if flip_mask==1
          surfMASK = surfMASK(end:-1:1,:,:);
       end
     end         

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

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

         s=sprintf('TEST = load_nii(''%s'');',[test_path source_name '_to_' target_name '.img']); 
         eval(s);
         TEST = double(TEST.img);
         %s=sprintf('fid = fopen(''%s'',''r'',''%s'');
         %',[test_path source_name '_to_' target_name '.img'],'ieee-be'); eval(s);
         %TEST = fread(fid,'int16');
         %s=sprintf('TEST = reshape(TEST,%s);',loaddims); eval(s);
         
       %-------------------------------------------------
       % Mask TEST and TRUE labels with intersecting TRUE 
       % (freesurfer) surface voxels
       %-------------------------------------------------
         if mask_surface == 1
               TEST = TEST .* surfMASK;
               TRUE = TRUE .* surfMASK;
       % max intersection with SPM_Normalize_direct_8mm 
       % (compensate for volume dimension difference) 
         elseif strcmp(reg_method,'SPM5_Normalize_direct_8mm')==1     
               %TEST = TEST(1:end-1,1:end-1,1:end-1);
               TEST = TEST(2:end,2:end,2:end); 
         end
                 
       %----------
       % Test plot
       %----------
         if plot_test==1
           ITRUE   = find(TRUE>0);
           [x,y,z] = ind2sub(size(TRUE),ITRUE);
           truexyz = [x,y,z];
           ITEST   = find(TEST>0);
           [x,y,z] = ind2sub(size(TEST),ITEST);
           testxyz = [x,y,z];
           for ii=1:size(truexyz,1)
               subplot(1,3,1), imagesc(TRUE(:,:,ii)); axis equal; axis off;
               subplot(1,3,2), imagesc(TEST(:,:,ii)); axis equal; axis off; 
               colormap(gray); 
               I=find(truexyz(:,3)==ii);
               subplot(1,3,3), plot(truexyz(I,2),truexyz(I,1),'ko'); hold on 
               I=find(testxyz(:,3)==ii);
                               plot(testxyz(I,2),testxyz(I,1),'r.'); hold off
                               axis([1 size(TRUE,1) 1 size(TRUE,2)]); axis equal; axis off;
               drawnow; pause(0);
           end
         end

       %-----------------------------------------
       % Evaluate labels assigned by the template
       %----------------------------------------- 
         [label_nTEST_nTRUE_nAND_nOR,...
          label_ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN,...
          nTEST_nTRUE_nAND_nOR,...
          ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN,...
          totalAND, totalOR, ndiff1] = evaluate_warps( TEST, TRUE, Labels, verbose );
      
          if save_results==1
            s='*nTEST_nTRUE* *ANDoverTRUE* totalAND totalOR ndiff1 Labels';  
            if mask_surface == 1
               eval(['save ' out_path source_name '_to_' target_name '_surf_eval.mat ' s ';']);
            else
               eval(['save ' out_path source_name '_to_' target_name '_eval.mat ' s ';']);
            end
          end
     end
   end
 end

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