wolffd@0: function [x0,y0] = centroid(x,y)
wolffd@0: % CENTROID Center of mass of a polygon.
wolffd@0: %	[X0,Y0] = CENTROID(X,Y) Calculates centroid 
wolffd@0: %	(center of mass) of planar polygon with vertices 
wolffd@0: %	coordinates X, Y.
wolffd@0: %	Z0 = CENTROID(X+i*Y) returns Z0=X0+i*Y0 the same
wolffd@0: %	as CENTROID(X,Y).
wolffd@0: 
wolffd@0: %  Copyright (c) 1995 by Kirill K. Pankratov,
wolffd@0: %       kirill@plume.mit.edu.
wolffd@0: %       06/01/95, 06/07/95
wolffd@0: 
wolffd@0: % Algorithm:
wolffd@0: %  X0 = Int{x*ds}/Int{ds}, where ds - area element
wolffd@0: %  so that Int{ds} is total area of a polygon.
wolffd@0: %    Using Green's theorem the area integral can be 
wolffd@0: %  reduced to a contour integral:
wolffd@0: %  Int{x*ds} = -Int{x^2*dy}, Int{ds} = Int{x*dy} along
wolffd@0: %  the perimeter of a polygon.
wolffd@0: %    For a polygon as a sequence of line segments
wolffd@0: %  this can be reduced exactly to a sum:
wolffd@0: %  Int{x^2*dy} = Sum{ (x_{i}^2+x_{i+1}^2+x_{i}*x_{i+1})*
wolffd@0: %  (y_{i+1}-y_{i})}/3;
wolffd@0: %  Int{x*dy} = Sum{(x_{i}+x_{i+1})(y_{i+1}-y_{i})}/2.
wolffd@0: %    Similarly
wolffd@0: %  Y0 = Int{y*ds}/Int{ds}, where
wolffd@0: %  Int{y*ds} = Int{y^2*dx} = 
wolffd@0: %  = Sum{ (y_{i}^2+y_{i+1}^2+y_{i}*y_{i+1})*
wolffd@0: %  (x_{i+1}-x_{i})}/3.
wolffd@0: 
wolffd@0:  % Handle input ......................
wolffd@0: if nargin==0, help centroid, return, end
wolffd@0: if nargin==1
wolffd@0:   sz = size(x);
wolffd@0:   if sz(1)==2      % Matrix 2 by n
wolffd@0:     y = x(2,:); x = x(1,:);
wolffd@0:   elseif sz(2)==2  % Matrix n by 2
wolffd@0:     y = x(:,2); x = x(:,1);
wolffd@0:   else
wolffd@0:     y = imag(x);
wolffd@0:     x = real(x);
wolffd@0:   end
wolffd@0: end 
wolffd@0: 
wolffd@0:  % Make a polygon closed ..............
wolffd@0: x = [x(:); x(1)];
wolffd@0: y = [y(:); y(1)];
wolffd@0: 
wolffd@0:  % Check length .......................
wolffd@0: l = length(x);
wolffd@0: if length(y)~=l
wolffd@0:   error(' Vectors x and y must have the same length')
wolffd@0: end
wolffd@0: 
wolffd@0:  % X-mean: Int{x^2*dy} ................
wolffd@0: del = y(2:l)-y(1:l-1);
wolffd@0: v = x(1:l-1).^2+x(2:l).^2+x(1:l-1).*x(2:l);
wolffd@0: x0 = v'*del;
wolffd@0: 
wolffd@0:  % Y-mean: Int{y^2*dx} ................
wolffd@0: del = x(2:l)-x(1:l-1);
wolffd@0: v = y(1:l-1).^2+y(2:l).^2+y(1:l-1).*y(2:l);
wolffd@0: y0 = v'*del;
wolffd@0: 
wolffd@0:  % Calculate area: Int{y*dx} ..........
wolffd@0: a = (y(1:l-1)+y(2:l))'*del;
wolffd@0: tol= 2*eps;
wolffd@0: if abs(a)<tol
wolffd@0:   disp(' Warning: area of polygon is close to 0')
wolffd@0:   a = a+sign(a)*tol+(~a)*tol;
wolffd@0: end
wolffd@0:  % Multiplier
wolffd@0: a = 1/3/a;
wolffd@0: 
wolffd@0:  % Divide by area .....................
wolffd@0: x0 = -x0*a;
wolffd@0: y0 =  y0*a;
wolffd@0: 
wolffd@0: if nargout < 2, x0 = x0+i*y0; end