polygonArea

 Compute the signed area of a polygon.

   A = polygonArea(POINTS);
   Compute area of a polygon defined by POINTS. POINTS is a N-by-2 array
   of double containing coordinates of vertices.
   
   Vertices of the polygon are supposed to be oriented Counter-Clockwise
   (CCW). In this case, the signed area is positive.
   If vertices are oriented Clockwise (CW), the signed area is negative.

   If polygon is self-crossing, the result is undefined.

   Examples
     % compute area of a simple shape
     poly = [10 10;30 10;30 20;10 20];
     area = polygonArea(poly)
     area = 
         200

     % compute area of CW polygon
     area2 = polygonArea(poly(end:-1:1, :))
     area2 = 
         -200

     % Computes area of a paper hen
     x = [0 10 20  0 -10 -20 -10 -10  0];
     y = [0  0 10 10  20  10  10  0 -10];
     poly = [x' y'];
     area = polygonArea(poly)
     area =
        400

     % Area of unit square with 25% hole
     pccw = [0 0; 1 0; 1 1; 0 1];
     pcw = pccw([1 4 3 2], :) * .5 + .25;
     polygonArea ([pccw; nan(1,2); pcw])
     ans =
        0.75

   References
   algo adapted from P. Bourke web page
   http://paulbourke.net/geometry/polygonmesh/

   See also:
   polygons2d, polygonCentroid, polygonSecondAreaMoments, triangleArea

0001 function area = polygonArea(poly, varargin)
0002 % Compute the signed area of a polygon.
0003 %
0004 %   A = polygonArea(POINTS);
0005 %   Compute area of a polygon defined by POINTS. POINTS is a N-by-2 array
0006 %   of double containing coordinates of vertices.
0007 %
0008 %   Vertices of the polygon are supposed to be oriented Counter-Clockwise
0009 %   (CCW). In this case, the signed area is positive.
0010 %   If vertices are oriented Clockwise (CW), the signed area is negative.
0011 %
0012 %   If polygon is self-crossing, the result is undefined.
0013 %
0014 %   Examples
0015 %     % compute area of a simple shape
0016 %     poly = [10 10;30 10;30 20;10 20];
0017 %     area = polygonArea(poly)
0018 %     area =
0019 %         200
0020 %
0021 %     % compute area of CW polygon
0022 %     area2 = polygonArea(poly(end:-1:1, :))
0023 %     area2 =
0024 %         -200
0025 %
0026 %     % Computes area of a paper hen
0027 %     x = [0 10 20  0 -10 -20 -10 -10  0];
0028 %     y = [0  0 10 10  20  10  10  0 -10];
0029 %     poly = [x' y'];
0030 %     area = polygonArea(poly)
0031 %     area =
0032 %        400
0033 %
0034 %     % Area of unit square with 25% hole
0035 %     pccw = [0 0; 1 0; 1 1; 0 1];
0036 %     pcw = pccw([1 4 3 2], :) * .5 + .25;
0037 %     polygonArea ([pccw; nan(1,2); pcw])
0038 %     ans =
0039 %        0.75
0040 %
0041 %   References
0042 %   algo adapted from P. Bourke web page
0043 %   http://paulbourke.net/geometry/polygonmesh/
0044 %
0045 %   See also:
0046 %   polygons2d, polygonCentroid, polygonSecondAreaMoments, triangleArea
0047 %
0048 
0049 %   ---------
0050 %   author : David Legland
0051 %   INRA - TPV URPOI - BIA IMASTE
0052 %   created the 05/05/2004.
0053 %
0054 
0055 %   HISTORY
0056 %   25/04/2005: add support for multiple polygons
0057 %   12/10/2007: update doc
0058 
0059 
0060 %% Process special cases
0061 
0062 % in case of polygon sets, computes the sum of polygon areas
0063 if iscell(poly)
0064     area = 0;
0065     for i = 1:length(poly)
0066         area = area + polygonArea(poly{i});
0067     end
0068     return;
0069 end
0070 
0071 % check there are enough points
0072 if size(poly, 1) < 2
0073     area = 0;
0074     return;
0075 end
0076 
0077 % case of polygons with holes -> computes the sum of areas
0078 if any(isnan(poly))
0079     area = sum(polygonArea(splitPolygons(poly)));
0080     return;
0081 end
0082 
0083 
0084 %% Process single polygons or single rings
0085 
0086 % extract coordinates
0087 if nargin == 1
0088     % polygon given as N-by-2 array
0089     px = poly(:, 1);
0090     py = poly(:, 2);
0091     
0092 elseif nargin == 2
0093     % poylgon given as two N-by-1 arrays
0094     px = poly;
0095     py = varargin{1};
0096 end
0097 
0098 % indices of next vertices
0099 N = length(px);
0100 iNext = [2:N 1];
0101 
0102 % compute area (vectorized version)
0103 area = sum(px .* py(iNext) - px(iNext) .* py) / 2;