function C = centroids(bnd)

%CENTROIDS takes a boundary file in standard format (below) and outputs 
%          the (approximate) centroid coordinates. Point coordinates are 
%          simply averaged, so that any straightline segments are discounted.

%Written by: TONY E. SMITH, April 21, 2002

% INPUT: X = (n:2) matrix describing boundaries of k polygons in the form:
%
%          [1   n1
%           x11 y11 
%           x12 y12
%           .
%           .
%           x11 y11
%           2   n2
%           x21 x22
%           .
%           .
%           .
%           k   nk
%           xk1 yk1
%           .
%           .
%           xk1 yk1]

%OUTPUTS:  C = (n:2) matrix of centroid coordinates


Z = polyform(bnd);

n = size(Z,1);

C = zeros(n:2);

for i = 1:n
   
   n1 = Z(i,1);
   
   n2 = Z(i,2);  %(n1,n2) define polygon i
   
   x_num = 0;
   y_num = 0;
   denom = 0;
   
   for j = n1:n2-1 %note: bnd(n2,:) = bnd(n1,:)
      
      x1 = bnd(j,1);
      y1 = bnd(j,2);
      x2 = bnd(j+1,1);
      y2 = bnd(j+1,2);
      
      len = sqrt((x1-x2)^2 + (y1-y2)^2);
      
      x_num = x_num + len*(x1 + x2)/2;
      y_num = y_num + len*(y1 + y2)/2;
      denom = denom + len;
      
   end
   
   C(i,1) = x_num/denom;
   C(i,2) = y_num/denom;   
      
end