Thanks to the many who responded to my question on how to calculate the area
inside a convex hull or polygon. Those responding were S.D. Byers, Barry
Rowlingson, Charles Berry, William Carlisle Thacker, Eric Turkheimer, Bill
Dunlap, Clint Bowman, Matt Kurbat, Wils Corrigan, Nick Ellis, Adrian
Baddeley, Duncan Mackay and Luciano, Molinari.
Bill Dunlap of Insightful Corporation provided S-Plus code that did
precisely what was needed. This code is quite similar to that provided by
Luciano Molinari, Nick Ellis and Barry Rowlingson. Responses are provided
below as 1) CODE 2) PSEUDO-CODE 3) REFERENCES and 4) SPATIAL MODULE.
Thanks again to all who responded. As always, I learned more and got more
from this list than I had expected.
Peter Scherer
Dow AgroSciences
Indianapolis
###
CODE
Bill Dunlap's function and explanation is presented here and is the code I
ended up using.
If x,y give the outline of a polygon then the following gives its
signed area (positive if the outline is clockwise, negative if
counterclockwise).
You can derive the formula by dividing the area into triangles, all
with a common vertex in an arbitrary position and using the cross product
formula for the area of a triangle based on the vectors for 2 sides.
You can also think of it as a discrete version of Green's theorem.
area <- function(x, y)
{
if(missing(y)) {
if(is.matrix(x) && ncol(x) == 2.) {
y <- x[, 2.]
x <- x[, 1.]
}
else if(!is.null(x$x) && !is.null(x$y)) {
y <- x$y
x <- x$x
}
}
x <- c(x, x[1.])
y <- c(y, y[1.])
i <- 2.:length(x)
return(0.5 * sum(x[i] * y[i - 1.] - x[i - 1.] * y[i]))
}
Nick Ellis provided a very nice set of S-Plus functions. Very similar to
that provided by Bill Dunlap.
Peter,
there was a question from Duncan Mackay on a related issue. I have attached
my reply to Duncan below for your information. That correspondence contains
some functions for calculating the area of a general polygon, which I
reproduce here for convenience.
area.in.polygon takes argument x which can be a list with components x and
y, a matrix with x in first comun and y in second or an object of class
"polygon" as created by as.polygon.
area.in.polygon <-
function(x)
{
x<-as.polygon(x)
n<-nrow(x)
D <- diff(x[c(1:n, 1), ])
A <- x[, 1] * D[, 2] - x[, 2] * D[, 1]
sum(A)/2
}
is.polygon <-
function(x)
{
!is.null(class(x)) && (class(x) == "polygon")
}
as.polygon <-
function(x)
{
if(is.polygon(x))
return(x)
if(is.matrix(x)) {
class(x) <- "polygon"
return(x)
}
x <- cbind(c(x$x, x$x[1]), c(x$y, x$y[1]))
class(x) <- "polygon"
attr(x, "closed") <- T
x
}
Nick Ellis
CSIRO Marine Research mailto:Nick.Ellis@marine.csiro.au
PO Box 120 ph +61 (07) 3826 7260
Cleveland QLD 4163 fax +61 (07) 3826 7222
Australia http://www.marine.csiro.au
Luciano Molinari provided a very similar function.
Below is a small function to this purpose.
No warranty whatsoever, but I hope it helps.
Luciano
------------------------
area
function(x, y)
{
if(is.matrix(x)) {
if(dim(x)[2] != 2)
stop("Error in input x")
y <- x[, 2]
x <- x[, 1]
}
x <- as.vector(x)
y <- as.vector(y)
if(length(x) != length(y))
stop("x and y must have the same length")
n <- length(x)
if(n <= 2)
0
else 0.5 * sum((x - c(x[3:n], x[1:2])) * c(y[2:n], y[1]))
}
Barry Rowlingson provided a general approach followed by some fortran code.
Calculate the sum of the areas defined by the trapeziums formed by pairs
of
points along the polygon and vertical lines to the X-axis. By
considering the
signs of the areas, the sum will add up to the area, since the area
below the
polygon is counted once positively and once negatively. Hmm, is that
clear?
Code perhaps... Some fortran?
function plarea(xp,yp,np)
c-------------------------------------------------------------------------
c
c find the area of the polygon defined by points in xp,yp
c
implicit real*8 (a-h,o-z)
dimension xp(np),yp(np)
totare=0
do is=1,np
x1=xp(is)
y1=yp(is)
if(is.eq.np)then
x2=xp(1)
y2=yp(1)
else
x2=xp(is+1)
y2=yp(is+1)
end if
c Find the area of the trapezium
totare = totare + (x2-x1)*(y2+y1)/2.0
end do
c return a positive value
plarea = abs(totare)
end
Easily incorporated into Splus, or just install my Splancs library and
you
have an 'area' function!
http://www.maths.lancs.ac.uk/~rowlings/Splancs
Barry
Charles Berry provided the following. I haven't yet further explored this
to see if chull can be used to calculate area. My examination of chull
indicated it could not, but perhaps I'm missing something.
chull(my.matrix), perhaps?
###
PSEUDO-CODE
Bill Thacker provided the following.
Isn't the formula for twice the area of a triangle with vertices a,b,
and c simply: xa(yb-yc) + xb(yc-ya) + xc(ya-yb) ?
If there is a point P within a polygon with vertices a, b, ..., you
can sum the areas of the triangles. The coordinates of P should drop
out, as it is chosen arbitrarily. For twice the area, I get:
xa*yb-xb*ya + xb*yc-xc*yb + ...
Clint Bowman's pseudo-code is as follows and provides a general sketch of
the method detailed elsewhere.
No code but it's pretty simple:
1) Pick any point within the convex hull (probably could be outside as well
but interior works better for the visualization.)
2) Each pair of consecutive points in the list that makes up the convex
hull,
when paired with the interior point, makes a triangle. You can easily
calculate the area knowing the coordinates of the three points.
3) Continue around the convex hull adding the triangular areas until you
reach
the starting point (I've forgotten whether the list has the starting point
repeated as the end point but a simple check will answer.)
4) Viola, you are done (Oh, wait, that's a musical joke - my apologies if
you
play that stringed instrument.)
Clint
Adrian Baddeley provided the following sketch.
The area inside an *arbitrary* closed polygon with vertices
v_1, ..., v_n with v_n = v_1 is the sum, over all edges e_i =
(v_i,v_{i+1}), of the signed area of Q(v_i, v_{i+1}), the
quadrangle enclosed by v_i, v_{i+1} and the projections of these
two points onto the x-axis. (The "discrete Stokes theorem")
It's easiest if the vertices are listed in anticlockwise
order around the polygon, and the x-axis is taken to be below
the entire polygon. Then the sign of the quadrangle area is
positive iff the x-coordinate of v_i is larger than the x-coordinate
of v_{i+1}.
###
REFERENCE
Eric Turkheimer provided a reference for some of his research.
I don't get to cite this very often... The first paper I ever wrote! It has
an algorithm....
Turkheimer, E., Yeo, R., & Bigler, E.D. (1983). Digital planimetry in
APLSF. Behavior Methods Research and Instrumentation, 15,
471-473.
Eric Turkheimer
Department of Psychology
Matt Kurbat provided the following reference.
Cormen, Leiserson, and Rivest's book (Introduction to Algorithms, MIT Press)
has pseudo-code, I believe, for area inside a convex hull.
Matt
Duncan Mackay provided the following link.
Here is a link to a nice computational geometry site that deal with polygon
calculations.
Duncan
http://compgeom.cs.uiuc.edu/~jeffe/compgeom/code.html#poly
###
SPATIAL MODULE
S.D Byers and Will Corrigan indicated that the Spatial module contained a
function for calculating the area of a polygon.
poly.area(x, y)
|