function [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 )          
%---------------------------------------------------------------------
% 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.
%
% Compare two sets of labels for an image volume.
% Determine the intersection and union label agreements 
% and type I and type II errors for template labels 
% assigned to independently labeled targets.
%
% Measure label agreements as intersection/true,
% intersection/union, and intersection/mean volume
% (the percent accord of Caviness et al (1996).
%
% INPUT:  TEST    = 3-D matrix: TEST (usually automated) labels
%         TRUE    = 3-D matrix: TRUE (manually assigned) labels
%         Labels  = list of labels (integers)
%         verbose = 1x1
%
% OUTPUT: label_nTEST_nTRUE_nAND_nOR:
%
%         [label, #TEST voxels, #TRUE voxels, #intersecting voxels, #union voxels]
%
%         label_ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN: 
%
%         [label, intersection/TRUE, intersection/union, intersection/mean volume,
%          false positive, false negative]
%
%         nTEST_nTRUE_nAND_nOR: 
%               same as label_... above summed over all brain regions
%         ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN: 
%               same as label... above calculated for the whole brain
%         totalAND = overall intersection of all Labels
%         totalOR  = overall union of all Labels
%         ndiff1   = number of voxels with a difference of one between
%                    TEST and TRUE labels
%
%--------------------------------------------------------------
% (c) 2008, @rno klein
%--------------------------------------------------------------

  format compact 
  plotem = 0;

%-------------------------------------------------------------
% Tally intersecting voxels with label values differing by one
%-------------------------------------------------------------
  I0diff = find( abs(TRUE-TEST)==1 );    
  ndiff1 = length(I0diff);

  Labels = unique(Labels(Labels>0));
  
%-------------------------------------------------
% Tally total intersection, union of TEST and TRUE
%-------------------------------------------------
  ITEST = find( ismember(TEST,Labels) );
  ITRUE = find( ismember(TRUE,Labels) );
  IAND  = intersect(ITEST, ITRUE);
  IOR   = union(ITEST, ITRUE);
  totalAND = size(IAND,1);
  totalOR  = size(IOR,1);

%-----------------------
% Tally totals per label
%----------------------- 
  label_nTEST_nTRUE_nAND_nOR                    = zeros(size(Labels,1),5);
  label_ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN = zeros(size(Labels,1),6);
  
  label_nTEST_nTRUE_nAND_nOR(:,1)                    = Labels;
  label_ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN(:,1) = Labels;

  for iLabel = 1 : length(Labels)

   ITRUE = find( TRUE == Labels(iLabel) );
   
   if size(ITRUE,1) > 0
        
        ITEST = find( TEST == Labels(iLabel) );

     %--------------------------------------------
     % Intersection, union of TEST and TRUE labels
     %--------------------------------------------
        IAND = intersect(ITEST, ITRUE);
        IOR  = union(ITEST, ITRUE);

        if plotem==1
           L = zeros(size(TRUE));
           L(ITRUE) = 1;
           L(ITEST) = L(ITEST)+1;
           islice = 125;
           imagesc(L(:,:,islice)); axis equal; axis off; drawnow; pause(1);
        end
   
        label_nTEST_nTRUE_nAND_nOR(iLabel,:) = [Labels(iLabel)...
                                   length(ITEST)...
                                   length(ITRUE)...
                                   length(IAND)...
                                   length(IOR)];
                          
        label_ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN(iLabel,:) = [Labels(iLabel)...
                                   length(IAND)/length(ITRUE)...
                                   length(IAND)/length(IOR)...
                                   length(IAND)*2/(length(ITEST)+length(ITRUE))...
                                   (length(IOR)-length(ITRUE))/length(ITEST)...
                                   (length(IOR)-length(ITEST))/length(ITRUE)];
     end
  end

  nTEST_nTRUE_nAND_nOR = sum(label_nTEST_nTRUE_nAND_nOR(:,2:5),1);

  ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN = [ nTEST_nTRUE_nAND_nOR(3)/nTEST_nTRUE_nAND_nOR(2)...
                                   nTEST_nTRUE_nAND_nOR(3)/nTEST_nTRUE_nAND_nOR(4)...
                                   nTEST_nTRUE_nAND_nOR(3)*2/(nTEST_nTRUE_nAND_nOR(1)+nTEST_nTRUE_nAND_nOR(2))...
                                   (nTEST_nTRUE_nAND_nOR(4)-nTEST_nTRUE_nAND_nOR(2))/nTEST_nTRUE_nAND_nOR(1)...
                                   (nTEST_nTRUE_nAND_nOR(4)-nTEST_nTRUE_nAND_nOR(1))/nTEST_nTRUE_nAND_nOR(2) ];

  if verbose > 1  
    %label_nTEST_nTRUE_nAND_nOR
    %label_ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN
    ANDoverTRUE_ANDoverOR_ANDovermean_FP_FN
    nTEST_nTRUE_nAND_nOR
    %totalAND
    %totalOR
    %ndiff1   
  end