Dear all,
is there any function to calculate the center of a polygon mass in R?
Actually I need to find the best location within polygons to place labels.
Thanks for any hint
Jens Oehlschl?gel
--
COMPUTERBILD 15/03: Premium-e-mail-Dienste im Test\ --------...{{dropped}}
Dear Jens,> -----Original Message----- > From: oehl_list at gmx.de [mailto:oehl_list at gmx.de] > Sent: 13 August 2003 10:01 > To: r-help at stat.math.ethz.ch > Subject: [R] placing labels in polygon center ? > > Dear all, > > is there any function to calculate the center of a polygon mass in R? > Actually I need to find the best location within polygons to > place labels.Trying to recall from my physics education, the mass centre is at the averages of the coordinates of the corners, isn't it? However, the task is more complex: your polygon might be X-shaped... HTH Thomas --- Thomas Hotz Research Associate in Medical Statistics University of Leicester United Kingdom Department of Epidemiology and Public Health 22-28 Princess Road West Leicester LE1 6TP Tel +44 116 252-5410 Fax +44 116 252-5423 Division of Medicine for the Elderly Department of Medicine The Glenfield Hospital Leicester LE3 9QP Tel +44 116 256-3643 Fax +44 116 232-2976
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.
I wrote:
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.
Here I was regarding a vector as a _sequence_.
The (one-step) rotation of c(u,v,w,x,y,z) is c(v,w,x,y,z,u).
This is pretty much the way APL uses the word "rotate" (the
vertical-bar-overstruck-with-a-circle operator).
Spencer Graves <spencer.graves at PDF.COM> replied:
I didn't study your code, but regarding a function to rotate a
vector: Multiplication by an orthogonal matrix does that.
This is a misunderstanding. We were both using the word "rotate" in
a standard way, the problem is that it has more than one "standard"
meaning.
As a matter of fact,
/ 0 1 0 0 \ / u \ / v \
| 0 0 1 0 | | v | | w |
| 0 0 0 1 | | w | = | x |
\ 1 0 0 0 / \ x / \ u /
so you *can* do the kind of rotation I want using a matrix multiplication,
and this is mathematically useful; it's just not a very good way to do
it in a computer.
This is rather simple-minded:
rot <- function(x, k=1) {
k <- k %% length(x)
x[c((k+1):length(x), 1:k)]
}
Andy
> -----Original Message-----
> From: Richard A. O'Keefe [mailto:ok at cs.otago.ac.nz]
> Sent: Thursday, August 14, 2003 1:28 AM
> To: r-help at stat.math.ethz.ch
> Subject: Re: [R] placing labels in polygon center ?
>
>
> I wrote:
> 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.
>
> Here I was regarding a vector as a _sequence_.
> The (one-step) rotation of c(u,v,w,x,y,z) is c(v,w,x,y,z,u).
> This is pretty much the way APL uses the word "rotate" (the
> vertical-bar-overstruck-with-a-circle operator).
>
> Spencer Graves <spencer.graves at PDF.COM> replied:
> I didn't study your code, but regarding a function to rotate a
> vector: Multiplication by an orthogonal matrix does that.
>
> This is a misunderstanding. We were both using the word
> "rotate" in a standard way, the problem is that it has more
> than one "standard" meaning.
>
> As a matter of fact,
> / 0 1 0 0 \ / u \ / v \
> | 0 0 1 0 | | v | | w |
> | 0 0 0 1 | | w | = | x |
> \ 1 0 0 0 / \ x / \ u /
> so you *can* do the kind of rotation I want using a matrix
> multiplication, and this is mathematically useful; it's just
> not a very good way to do it in a computer.
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo> /r-help
>
------------------------------------------------------------------------------
Notice: This e-mail message, together with any attachments, contains
information of Merck & Co., Inc. (Whitehouse Station, New Jersey, USA),
and/or
its affiliates (which may be known outside the United States as Merck Frosst,
Merck Sharp & Dohme or MSD) that may be confidential, proprietary
copyrighted
and/or legally privileged, and is intended solely for the use of the
individual or entity named on this message. If you are not the intended
recipient, and have received this message in error, please immediately return
this by e-mail and then delete it.
"Petr Pikal" <petr.pikal at precheza.cz> wrote (embedded in much
XML):
Not sure about efficency and it is not very general solution but
you maybe can use embed() function
It's an interesting suggestion, but I don't see *how*.
Here's the function I want, only I'd like it to be something built in.
I've limited it to the where the rotation count is non-negative, just
to keep it simple. For the same reason, I've limited it to vectors,
although the operation can be generalised in useful ways to matrices
and arrays with any number of subscripts.
rotate <- function (x, count=1) {
x[((0:(length(x)-1)) + (count %% length(x))) %% length(x) + 1]
}
x <- 1:5
rotate(x)
=> 2 3 4 5 1
rotate(x, 2)
=> 3 4 5 1 2
rotate(x, 3)
=> 4 5 1 2 3
One important thing about this is that the result of rotate(x)
is a permutation of x; it has the same number of elements, the same
elements, only the order is changed.
The reason I don't see _how_ to use embed() to do this is that the
result of embed is a matrix (but I want a vector) which has fewer rows
that the original (but I want the same number of elements).
embed(x, n) => a matrix with length(x)-n+1 rows and n columns.
In the polygon application, embed(xs, 2) is very close to what's needed.
If the input to the polygon functions had the first element repeated at
the end, it would be just right:
a
b => b a
c c b
a a c
The problem is that the input isn't a b c a, it's a b c.
Let's not take up too many people's time with this, OK?
I hoped there might be a built-in function I had missed; apparently
there isn't. End of story.
Barry Rowlingson <B.Rowlingson at lancaster.ac.uk> wrote:
Do you want:
c(x[-1],x[1]) for a one-step 'rotation'?
That's just the kind of thing I did, except that it's ugly.
I've browsed src/main/subscript.c and got rather lost, but
it looks very much as though x[-1] starts by making a
logical vector c(FALSE,TRUE,...,TRUE) and then that's used
as an index to make a copy of (part of) x, then there's
x[1] (in general, x[1:n]) to be copied, and finally these
copies are pasted together. That's 2 or 3 times as much
memory allocated as is actually kept; and it's pretty obvious
that it can be done in one quick pass with no redundant memory
allocation at the implementation level.
With rotation being so useful in APL, I rather expected that
there would be some C level implementation of rotate() or of
something better. Note that rotate() should really have
three arguments:
rotate(array, amount=1, axis=1)