From octave-sources-request at bevo dot che dot wisc dot edu  Wed Apr 14 15:05:06 1999
Subject: comp geom
From: "John W. Eaton" <jwe at bevo dot che dot wisc dot edu>
To: doolin at alpha-8 dot CE dot Berkeley dot EDU
Cc: octave-sources at bevo dot che dot wisc dot edu
Date: Wed, 14 Apr 1999 15:04:46 -0500 (CDT)

On 14-Apr-1999, doolin at alpha-8 dot CE dot Berkeley dot EDU <doolin@alpha-8.CE.Berkeley.EDU> wrote:

| Another addition to a computation geometry 
| tool set.
| 
| As usual, suggestions, comments etc encouraged.


| 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

Here is another way to write this that avoids the loop, which may make
it faster if the number of rows of p is large.

  function a = polyarea (p)
    if (nargin == 1)
      [r, c] = size (p);
      if (c == 2 && r > 3)
	a = sum (sum (p .* [p(2:r,2), -p(2:r,1); p(1,2), -p(1,1)])) / 2;
      else
        error ("polyarea: expecting matrix with 2 columns and 3 or more rows");
      endif
    else
      usage ("polyarea (p)");
    endif
  endfunction

In addition to checking to see that the segments don't cross (as you
mention in your comments) it would probably also be good (but harder)
to check to see that the vertices are given in the correct order since it
seems that they need to be specified in a counter-clockwise order.
For example,

  polyarea ([1,1; -1,1; -1,-1; 1,-1])  ==>   4
  polyarea ([1,1; 1,-1; -1,-1; -1,1])  ==>  -4

jwe