From octave-sources-request at bevo dot che dot wisc dot edu  Wed Apr 14 10:53:38 1999
Subject: comp geom
From: doolin at alpha-8 dot CE dot Berkeley dot EDU
To: octave-sources at bevo dot che dot wisc dot edu
Cc: doolin at alpha-8 dot CE dot Berkeley dot EDU
Date: Wed, 14 Apr 1999 08:53:35 -0700


Another addition to a computation geometry 
tool set.

As usual, suggestions, comments etc encouraged.



% [ area ] = polyarea(polygon)
%
% Determines area of a polygon by triangle method.
%
% input:
%   polygon:  Matrix of x (1st column), y (2d column) vertex
%             pairs (rows).
%
% output:
%   area:  Area of polygon.  Should be positive.
%
% todo:  Add moments for centroid, etc.
%
% bugs and limitations:  
%        Probably ought to be an optional check to make sure that
%        traversing the vertices doesn't make any sides cross 
%        (Is simple closed curve the technical definition of this?). 
%        
% This script properly belongs in octave/script/compgeom/
% (computational geometry).
%
% This code is copyright David M. Doolin under the GNU
% General Public License (GPL). This program is distributed 
% in the hope that it will be useful, but WITHOUT ANY WARRANTY; 
% without even the implied warranty of MERCHANTABILITY or 
% FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public 
% License for more details.  You should have received a 
% copy of the GNU General Public License along with this file.  
% If not, write to the Free Software Foundation,
% 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
%
% $Author: doolin $  doolin at ce dot berkeley dot edu
% $Date: 1999/04/14 15:52:20 $
% $Source: /shag/homes/doolin/cvsroot/octave/compgeom/polyarea.m,v $
% $Revision: 1.1 $

function [ area ] = polyarea(p)

if (nargin != 1)
   usage("polyarea(polygon)");
end

%  For matlab compatibility
[numvertices scratch] = size(p);

det = 0;

% Save the initial vertex for the last triangle
% f == firstrow
f = p(1,:);

for j = 1:numvertices-1
   det = det + p(j,1)*p(j+1,2) - p(j+1,1)*p(j,2);
end
% We still need j.
j = j + 1;
det = det +  p(j,1)*f(1,2) - f(1,1)*p(j,2);

area = det/2.0;

endfunction