|
wolffd@0
|
1 function [x0,y0] = centroid(x,y)
|
|
wolffd@0
|
2 % CENTROID Center of mass of a polygon.
|
|
wolffd@0
|
3 % [X0,Y0] = CENTROID(X,Y) Calculates centroid
|
|
wolffd@0
|
4 % (center of mass) of planar polygon with vertices
|
|
wolffd@0
|
5 % coordinates X, Y.
|
|
wolffd@0
|
6 % Z0 = CENTROID(X+i*Y) returns Z0=X0+i*Y0 the same
|
|
wolffd@0
|
7 % as CENTROID(X,Y).
|
|
wolffd@0
|
8
|
|
wolffd@0
|
9 % Copyright (c) 1995 by Kirill K. Pankratov,
|
|
wolffd@0
|
10 % kirill@plume.mit.edu.
|
|
wolffd@0
|
11 % 06/01/95, 06/07/95
|
|
wolffd@0
|
12
|
|
wolffd@0
|
13 % Algorithm:
|
|
wolffd@0
|
14 % X0 = Int{x*ds}/Int{ds}, where ds - area element
|
|
wolffd@0
|
15 % so that Int{ds} is total area of a polygon.
|
|
wolffd@0
|
16 % Using Green's theorem the area integral can be
|
|
wolffd@0
|
17 % reduced to a contour integral:
|
|
wolffd@0
|
18 % Int{x*ds} = -Int{x^2*dy}, Int{ds} = Int{x*dy} along
|
|
wolffd@0
|
19 % the perimeter of a polygon.
|
|
wolffd@0
|
20 % For a polygon as a sequence of line segments
|
|
wolffd@0
|
21 % this can be reduced exactly to a sum:
|
|
wolffd@0
|
22 % Int{x^2*dy} = Sum{ (x_{i}^2+x_{i+1}^2+x_{i}*x_{i+1})*
|
|
wolffd@0
|
23 % (y_{i+1}-y_{i})}/3;
|
|
wolffd@0
|
24 % Int{x*dy} = Sum{(x_{i}+x_{i+1})(y_{i+1}-y_{i})}/2.
|
|
wolffd@0
|
25 % Similarly
|
|
wolffd@0
|
26 % Y0 = Int{y*ds}/Int{ds}, where
|
|
wolffd@0
|
27 % Int{y*ds} = Int{y^2*dx} =
|
|
wolffd@0
|
28 % = Sum{ (y_{i}^2+y_{i+1}^2+y_{i}*y_{i+1})*
|
|
wolffd@0
|
29 % (x_{i+1}-x_{i})}/3.
|
|
wolffd@0
|
30
|
|
wolffd@0
|
31 % Handle input ......................
|
|
wolffd@0
|
32 if nargin==0, help centroid, return, end
|
|
wolffd@0
|
33 if nargin==1
|
|
wolffd@0
|
34 sz = size(x);
|
|
wolffd@0
|
35 if sz(1)==2 % Matrix 2 by n
|
|
wolffd@0
|
36 y = x(2,:); x = x(1,:);
|
|
wolffd@0
|
37 elseif sz(2)==2 % Matrix n by 2
|
|
wolffd@0
|
38 y = x(:,2); x = x(:,1);
|
|
wolffd@0
|
39 else
|
|
wolffd@0
|
40 y = imag(x);
|
|
wolffd@0
|
41 x = real(x);
|
|
wolffd@0
|
42 end
|
|
wolffd@0
|
43 end
|
|
wolffd@0
|
44
|
|
wolffd@0
|
45 % Make a polygon closed ..............
|
|
wolffd@0
|
46 x = [x(:); x(1)];
|
|
wolffd@0
|
47 y = [y(:); y(1)];
|
|
wolffd@0
|
48
|
|
wolffd@0
|
49 % Check length .......................
|
|
wolffd@0
|
50 l = length(x);
|
|
wolffd@0
|
51 if length(y)~=l
|
|
wolffd@0
|
52 error(' Vectors x and y must have the same length')
|
|
wolffd@0
|
53 end
|
|
wolffd@0
|
54
|
|
wolffd@0
|
55 % X-mean: Int{x^2*dy} ................
|
|
wolffd@0
|
56 del = y(2:l)-y(1:l-1);
|
|
wolffd@0
|
57 v = x(1:l-1).^2+x(2:l).^2+x(1:l-1).*x(2:l);
|
|
wolffd@0
|
58 x0 = v'*del;
|
|
wolffd@0
|
59
|
|
wolffd@0
|
60 % Y-mean: Int{y^2*dx} ................
|
|
wolffd@0
|
61 del = x(2:l)-x(1:l-1);
|
|
wolffd@0
|
62 v = y(1:l-1).^2+y(2:l).^2+y(1:l-1).*y(2:l);
|
|
wolffd@0
|
63 y0 = v'*del;
|
|
wolffd@0
|
64
|
|
wolffd@0
|
65 % Calculate area: Int{y*dx} ..........
|
|
wolffd@0
|
66 a = (y(1:l-1)+y(2:l))'*del;
|
|
wolffd@0
|
67 tol= 2*eps;
|
|
wolffd@0
|
68 if abs(a)<tol
|
|
wolffd@0
|
69 disp(' Warning: area of polygon is close to 0')
|
|
wolffd@0
|
70 a = a+sign(a)*tol+(~a)*tol;
|
|
wolffd@0
|
71 end
|
|
wolffd@0
|
72 % Multiplier
|
|
wolffd@0
|
73 a = 1/3/a;
|
|
wolffd@0
|
74
|
|
wolffd@0
|
75 % Divide by area .....................
|
|
wolffd@0
|
76 x0 = -x0*a;
|
|
wolffd@0
|
77 y0 = y0*a;
|
|
wolffd@0
|
78
|
|
wolffd@0
|
79 if nargout < 2, x0 = x0+i*y0; end
|
