[R] placing labels in polygon center ?

[Richard A. O'Keefe]

Barry Rowlingson <B.Rowlingson at lancaster.ac.uk> provided
functions PolygonArea and PolygonCenterOfMass.

As an exercise in R programming, I thought "why don't I vectorise these
and then see if it makes a practical difference".

Here are my versions of his functions.  Somehow I ended up with a sign
error when I entered his centroid code, so I had better exhibit the code
that I actually tested.

polygon.area <- function (polygon) {
    N <- dim(polygon)[1]
    area <- 0
    for (i in 1:N) {
       j <- i %% N + 1
       area <- area + polygon[i,1]*polygon[j,2] - polygon[i,2]*polygon[j,1]
   }
   abs(area/2)
}

polygon.centroid <- function(polygon) {
    N <- dim(polygon)[1]
    cx <- cy <- 0
    for (i in 1:N) {
        j <- i %% N + 1
        factor <- polygon[j,1]*polygon[i,2] - polygon[i,1]*polygon[j,2]
        cx <- cx + (polygon[i,1]+polygon[j,1])*factor
        cy <- cy + (polygon[i,2]+polygon[j,2])*factor
    }
    factor <- 1/(6*polygon.area(polygon))
    c(cx*factor, cy*factor)
}
 
Here are vectorised versions.  I found myself wishing for a function to
rotate a vector.  Is there one?  I know about ?lag, but
help.search("rotate") didn't find anything to the point.

vectorised.area <- function(polygon) {
    ix <- c(2:dim(polygon)[1], 1)
    xs <- polygon[,1] 
    ys <- polygon[,2]
    abs(sum(xs*ys[ix]) - sum(xs[ix]*ys))/2
}
 
vectorised.centroid <- function(polygon) {
    ix <- c(2:dim(polygon)[1], 1)
    xs <- polygon[,1]; xr <- xs[ix]
    ys <- polygon[,2]; yr <- ys[ix]
    factor <- xr*ys - xs*yr
    cx <- sum((xs+xr)*factor)
    cy <- sum((ys+yr)*factor)
    scale <- 3*abs(sum(xs*yr) - sum(xr*ys))
    c(cx/scale, cy/scale)
}

Test case 1: unit square.

> p <- rbind(c(0,0), c(0,1), c(1,1), c(1,0))
> polygon.area(p)
[1] 1
> vectorised.area(p)
[1] 1
> polygon.centroid(p)
[1] 0.5 0.5
> vectorised.centroid(p)
[1] 0.5 0.5
> system.time(for (i in 1:1000) polygon.area(p))
[1] 0.56 0.02 0.58 0.00 0.00
> system.time(for (i in 1:1000) vectorised.area(p))
[1] 0.22 0.03 0.25 0.00 0.00
> system.time(for (i in 1:1000) polygon.centroid(p))
[1] 1.56 0.06 1.66 0.00 0.00
> system.time(for (i in 1:1000) vectorised.centroid(p))
[1] 0.35 0.04 0.39 0.00 0.00

Even for a polygon this small, vectorising pays off.

Test case 2: random 20-gon.

> p <- cbind(runif(20), runif(20)) 
> polygon.area(p)
[1] 0.2263327
> vectorised.area(p)
[1] 0.2263327
> polygon.centroid(p)
[1] 0.6820708 0.5196700
> vectorised.centroid(p)
[1] 0.6820708 0.5196700
> system.time(for (i in 1:1000) polygon.area(p))
[1] 2.49 0.03 2.61 0.00 0.00
> system.time(for (i in 1:1000) vectorised.area(p))
[1] 0.29 0.05 0.34 0.00 0.00
> system.time(for (i in 1:1000) polygon.centroid(p))
[1] 7.29 0.07 7.70 0.00 0.00
> system.time(for (i in 1:1000) vectorised.centroid(p))
[1] 0.45 0.05 0.51 0.00 0.00

I was expecting the 20-gon version to be faster; what I did not expect
was that vectorising would pay off even for a quadrilateral.  In fact,

> p <- rbind(c(0,0), c(0,1), c(1,0))
> system.time(for (i in 1:1000) polygon.centroid(p))
[1] 1.25 0.04 1.31 0.00 0.00
> system.time(for (i in 1:1000) vectorised.centroid(p))
[1] 0.33 0.07 0.40 0.00 0.00

it even pays off for a triangle.