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