function [Boundaries,Borders] = make_boundaries( LABELS, Labels, dimensions, verbose );
%----------------------------------------------------------------------
% This program was used to generate the label boundaries
% 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.
%
% Create a 3-D boundary for each labeled region of the template,
%
% NOTE: Coregister the template labels with the template image.
%
% INPUT:  LABELS = 3-D matrix of labeled voxels (template dimensions)
%                      
% OUTPUT: Boundaries: Nx4 matrix: [x,y,z,label]
%         Borders:    Nx4 matrix: [x,y,z,label]
%
% (c) 2008, @rno klein
%--------------------------------------------------------------------------

  neighborhood_distance = 2; % Integer above maximum distance defining boundary
  plot_test             = 0;
  boundaries            = 0;
  
%--------------------------------------
% Determine the boundary for each label
%--------------------------------------
  if boundaries==1
       BOUNDARIES = zeros(dimensions);
  else Boundaries = [];
  end
  BORDERS    = zeros(dimensions);

  for iL = 1 : length(Labels)

     Label = Labels(iL);

     IL = find( LABELS == Label );

     if verbose > 0
        disp(['    Constructing boundary for label ' num2str(Label)]);
     end

     if length(IL) > 0

      % full_surface
        if boundaries==1
           D = bwdist( (LABELS~=Label), 'cityblock' );
           I = find( (D>0) & (D<neighborhood_distance) & (LABELS==Label) );

           BOUNDARIES(I) = Label;
        end
      
      % borders
        D = bwdist( (LABELS~=Label & LABELS>0), 'cityblock' );
        I = find( (D>0) & (D<neighborhood_distance) & (LABELS==Label) );

        BORDERS(I) = Label;

     end   

  end

  if boundaries==1
     IBOU             = find( BORDERS == 0 & BOUNDARIES > 0 );
  end
  IBOR             = find( BORDERS  > 0 );
  
  if boundaries==1
     [xbou ybou zbou] = ind2sub( dimensions, IBOU );
  end
  [xbor ybor zbor] = ind2sub( dimensions, IBOR );

  if boundaries==1
     Boundaries       = [xbou ybou zbou BOUNDARIES(IBOU)];
  end
  Borders          = [xbor ybor zbor BORDERS(IBOR)];

  if plot_test == 1
     figure
     if boundaries==1
        plot3(xbou,ybou,zbou,'b.'); hold on
     end
     plot3(xbor,ybor,zbor,'r.'); axis equal;
     title('Template boundaries (blue) and borders (red)');
     drawnow
  end