function [cx,cy]=PolygonCentroid(p)

%[cx,cy]=PolygonCentroid(p) returns the centroid of a polygon

%cx=\frac{1}{6A}/sum_0^{n-1} (x_i+x_{i+1})(x_i y_{i+1}-x_{i+1}y_i)
%cy=\frac{1}{6A}/sum_0^{n-1} (y_i+y_{i+1})(x_i y_{i+1}-x_{i+1}y_i)


A=(SignedArea(p));

cx=1/(6*A)*...
    sum(...
    (p(1,:)+[p(1,2:end),p(1,1)]).*...
    (...
    p(1,:).*[p(2,2:end),p(2,1)]-[p(1,2:end),p(1,1)].*p(2,:)...
    )...
    );

cy=1/(6*A)*...
    sum(...
    (p(2,:)+[p(2,2:end),p(2,1)]).*...
    (...
    p(1,:).*[p(2,2:end),p(2,1)]-[p(1,2:end),p(1,1)].*p(2,:)...
    )...
    );


if nargout < 2
    cx = [cx;cy];
end