Re: [R] Faster union of polygons?
Martin,
thanks a lot! This speeds things up so much
Do you mind if I bundle your punion function in a package I am
developing (but of course I will name you the author of the function)?
greetings,
Remko
On Fri, Jun 4, 2010 at 12:17 PM, Martin Morgan wrote:
> On 06/03/2010 04:54 PM, Remko Duursma wrote:
>> Thanks for the tip - this cleans up the code a lot!
>> Unfortunately, there is no gain in speed.
>
> Playing a little bit dirty,
>
> punion <-
> function(...)
> {
> n <- nargs()
> if (0L == n) new("gpc.poly")
> else if (1L == n && is(..1, "gpc.poly")) ..1
> else {
> polygons <- list(...)
> if (!all(sapply(polygons, is, "gpc.poly")))
> stop("'...' must all be 'gpc.poly'")
> ## avoid method look-up
> to_numeric <- selectMethod(coerce, c("gpc.poly", "numeric"))
> vec <- to_numeric(polygons[[1]])
> for (p in polygons[-1]) {
> clip <- to_numeric(p)
> vec <- .Call("Rgpc_polygon_clip", vec, clip, 3,
> PACKAGE="gpclib")
> }
> if (identical(vec, 0))
> new("gpc.poly")
> else
> as(vec, "gpc.poly")
> }
> }
>
> is about 4x faster on your example
>
>> replicate(5, system.time(Reduce(union, leaves)))
> [,1] [,2] [,3] [,4] [,5]
> user.self 1.272 1.272 1.272 1.268 1.268
> sys.self 0.000 0.000 0.000 0.000 0.000
> elapsed 1.271 1.272 1.272 1.273 1.281
> user.child 0.000 0.000 0.000 0.000 0.000
> sys.child 0.000 0.000 0.000 0.000 0.000
>> replicate(5, system.time(do.call(punion, leaves)))
> [,1] [,2] [,3] [,4] [,5]
> user.self 0.308 0.312 0.304 0.308 0.312
> sys.self 0.004 0.000 0.004 0.004 0.000
> elapsed 0.311 0.311 0.309 0.314 0.317
> user.child 0.000 0.000 0.000 0.000 0.000
> sys.child 0.000 0.000 0.000 0.000 0.000
>
> Rprof suggests that most of the time is now in the C code
>
>> Rprof("/tmp/leaves.Rprof")
>> x <- replicate(5, system.time(do.call(punion, leaves)))
>> Rprof(NULL)
>> summaryRprof("/tmp/leaves.Rprof")
> $by.self
> self.time self.pct total.time total.pct
> ".Call" 1.24 69.7 1.24 69.7
> "gc" 0.24 13.5 0.24 13.5
> "FUN" 0.08 4.5 1.78 100.0
> [...SNIP...]
>
> Martin
>
>>
>> remko
>>
>>
>>
>> On Thu, Jun 3, 2010 at 10:46 PM, nikhil kaza wrote:
>>> Reduce might work. Not sure about the speed advantages though. It does
>>> simplify code.
>>>
>>> Unionall <- function(x) Reduce('union', x)
>>> leaveout <- Unionall(leaves)
>>>
>>>
>>> On Tue, Jun 1, 2010 at 9:53 PM, Remko Duursma
>>> wrote:
Dear R-helpers,
thanks for yesterday's speeding-up tip. Here is my next query:
I have lots of polygons (not necessarily convex ones, and they never
have holes) given by x,y coordinates.
I want to get the polygon that is the union of these polygons. This is
my current method, but I am hoping there is a faster method (up to
thousands of polygons, each with ca. 40 xy points).
Example:
library(gpclib)
# A polygon
leaf <- structure(c(0, 1, 12.9, 16.5, 18.8, 17, 16.8, 15.5, 12.1, 8.2,
6.3, 5, 2, 0, -1.5, -4.3, -6.6, -10.3, -14.8, -19.4, -22.2, -23.5,
-22.2, -17.6, -7.8, 0, 0, -2.4, 2.8, 8.9, 19.9, 33.9, 34.8, 40.4,
49.7, 69.2, 77.4, 83.4, 91.4, 99, 92.8, 87.3, 81.2, 71.1, 57.6,
45.4, 39.2, 26, 15.6, 5.3, 0.6, 0), .Dim = c(26L, 2L), .Dimnames = list(
NULL, c("X", "Y")))
# Lots of polygons:
releaf <-
function(leaf)cbind(leaf[,1]+rnorm(1,0,50),leaf[,2]+rnorm(1,0,50))
leaves <- replicate(500, releaf(leaf), simplify=FALSE)
# Make into gpc.poly class:
leaves <- lapply(leaves, as, "gpc.poly")
# Make union .
system.time({
leavesoutline <- union(leaves[[1]], leaves[[2]])
for(i in 3:length(leaves))leavesoutline <- union(leavesoutline,
leaves[[i]])
})
# about 1sec here.
# Check it:
plot(leavesoutline)
thanks!
Remko
-
Remko Duursma
Research Lecturer
Centre for Plants and the Environment
University of Western Sydney
Hawkesbury Campus
Richmond NSW 2753
Dept of Biological Science
Macquarie University
North Ryde NSW 2109
Australia
Mobile: +61 (0)422 096908
www.remkoduursma.com
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
>>>
>>>
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read
Re: [R] Faster union of polygons?
On 06/03/2010 04:54 PM, Remko Duursma wrote:
> Thanks for the tip - this cleans up the code a lot!
> Unfortunately, there is no gain in speed.
Playing a little bit dirty,
punion <-
function(...)
{
n <- nargs()
if (0L == n) new("gpc.poly")
else if (1L == n && is(..1, "gpc.poly")) ..1
else {
polygons <- list(...)
if (!all(sapply(polygons, is, "gpc.poly")))
stop("'...' must all be 'gpc.poly'")
## avoid method look-up
to_numeric <- selectMethod(coerce, c("gpc.poly", "numeric"))
vec <- to_numeric(polygons[[1]])
for (p in polygons[-1]) {
clip <- to_numeric(p)
vec <- .Call("Rgpc_polygon_clip", vec, clip, 3,
PACKAGE="gpclib")
}
if (identical(vec, 0))
new("gpc.poly")
else
as(vec, "gpc.poly")
}
}
is about 4x faster on your example
> replicate(5, system.time(Reduce(union, leaves)))
[,1] [,2] [,3] [,4] [,5]
user.self 1.272 1.272 1.272 1.268 1.268
sys.self 0.000 0.000 0.000 0.000 0.000
elapsed1.271 1.272 1.272 1.273 1.281
user.child 0.000 0.000 0.000 0.000 0.000
sys.child 0.000 0.000 0.000 0.000 0.000
> replicate(5, system.time(do.call(punion, leaves)))
[,1] [,2] [,3] [,4] [,5]
user.self 0.308 0.312 0.304 0.308 0.312
sys.self 0.004 0.000 0.004 0.004 0.000
elapsed0.311 0.311 0.309 0.314 0.317
user.child 0.000 0.000 0.000 0.000 0.000
sys.child 0.000 0.000 0.000 0.000 0.000
Rprof suggests that most of the time is now in the C code
> Rprof("/tmp/leaves.Rprof")
> x <- replicate(5, system.time(do.call(punion, leaves)))
> Rprof(NULL)
> summaryRprof("/tmp/leaves.Rprof")
$by.self
self.time self.pct total.time total.pct
".Call" 1.24 69.7 1.24 69.7
"gc" 0.24 13.5 0.24 13.5
"FUN" 0.08 4.5 1.78 100.0
[...SNIP...]
Martin
>
> remko
>
>
>
> On Thu, Jun 3, 2010 at 10:46 PM, nikhil kaza wrote:
>> Reduce might work. Not sure about the speed advantages though. It does
>> simplify code.
>>
>> Unionall <- function(x) Reduce('union', x)
>> leaveout <- Unionall(leaves)
>>
>>
>> On Tue, Jun 1, 2010 at 9:53 PM, Remko Duursma
>> wrote:
>>>
>>> Dear R-helpers,
>>>
>>> thanks for yesterday's speeding-up tip. Here is my next query:
>>>
>>> I have lots of polygons (not necessarily convex ones, and they never
>>> have holes) given by x,y coordinates.
>>>
>>> I want to get the polygon that is the union of these polygons. This is
>>> my current method, but I am hoping there is a faster method (up to
>>> thousands of polygons, each with ca. 40 xy points).
>>>
>>> Example:
>>>
>>> library(gpclib)
>>>
>>> # A polygon
>>> leaf <- structure(c(0, 1, 12.9, 16.5, 18.8, 17, 16.8, 15.5, 12.1, 8.2,
>>> 6.3, 5, 2, 0, -1.5, -4.3, -6.6, -10.3, -14.8, -19.4, -22.2, -23.5,
>>> -22.2, -17.6, -7.8, 0, 0, -2.4, 2.8, 8.9, 19.9, 33.9, 34.8, 40.4,
>>> 49.7, 69.2, 77.4, 83.4, 91.4, 99, 92.8, 87.3, 81.2, 71.1, 57.6,
>>> 45.4, 39.2, 26, 15.6, 5.3, 0.6, 0), .Dim = c(26L, 2L), .Dimnames = list(
>>>NULL, c("X", "Y")))
>>>
>>> # Lots of polygons:
>>> releaf <-
>>> function(leaf)cbind(leaf[,1]+rnorm(1,0,50),leaf[,2]+rnorm(1,0,50))
>>> leaves <- replicate(500, releaf(leaf), simplify=FALSE)
>>>
>>> # Make into gpc.poly class:
>>> leaves <- lapply(leaves, as, "gpc.poly")
>>>
>>> # Make union .
>>> system.time({
>>> leavesoutline <- union(leaves[[1]], leaves[[2]])
>>> for(i in 3:length(leaves))leavesoutline <- union(leavesoutline,
>>> leaves[[i]])
>>> })
>>> # about 1sec here.
>>>
>>> # Check it:
>>> plot(leavesoutline)
>>>
>>>
>>>
>>> thanks!
>>>
>>> Remko
>>>
>>>
>>> -
>>> Remko Duursma
>>> Research Lecturer
>>>
>>> Centre for Plants and the Environment
>>> University of Western Sydney
>>> Hawkesbury Campus
>>> Richmond NSW 2753
>>>
>>> Dept of Biological Science
>>> Macquarie University
>>> North Ryde NSW 2109
>>> Australia
>>>
>>> Mobile: +61 (0)422 096908
>>> www.remkoduursma.com
>>>
>>> __
>>> R-help@r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>
>>
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
Martin Morgan
Computational Biology / Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N.
PO Box 19024 Seattle, WA 98109
Location: Arnold Building M1 B861
Phone: (206) 667-2793
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r
Re: [R] Faster union of polygons?
Thanks for the tip - this cleans up the code a lot!
Unfortunately, there is no gain in speed.
remko
On Thu, Jun 3, 2010 at 10:46 PM, nikhil kaza wrote:
> Reduce might work. Not sure about the speed advantages though. It does
> simplify code.
>
> Unionall <- function(x) Reduce('union', x)
> leaveout <- Unionall(leaves)
>
>
> On Tue, Jun 1, 2010 at 9:53 PM, Remko Duursma
> wrote:
>>
>> Dear R-helpers,
>>
>> thanks for yesterday's speeding-up tip. Here is my next query:
>>
>> I have lots of polygons (not necessarily convex ones, and they never
>> have holes) given by x,y coordinates.
>>
>> I want to get the polygon that is the union of these polygons. This is
>> my current method, but I am hoping there is a faster method (up to
>> thousands of polygons, each with ca. 40 xy points).
>>
>> Example:
>>
>> library(gpclib)
>>
>> # A polygon
>> leaf <- structure(c(0, 1, 12.9, 16.5, 18.8, 17, 16.8, 15.5, 12.1, 8.2,
>> 6.3, 5, 2, 0, -1.5, -4.3, -6.6, -10.3, -14.8, -19.4, -22.2, -23.5,
>> -22.2, -17.6, -7.8, 0, 0, -2.4, 2.8, 8.9, 19.9, 33.9, 34.8, 40.4,
>> 49.7, 69.2, 77.4, 83.4, 91.4, 99, 92.8, 87.3, 81.2, 71.1, 57.6,
>> 45.4, 39.2, 26, 15.6, 5.3, 0.6, 0), .Dim = c(26L, 2L), .Dimnames = list(
>> NULL, c("X", "Y")))
>>
>> # Lots of polygons:
>> releaf <-
>> function(leaf)cbind(leaf[,1]+rnorm(1,0,50),leaf[,2]+rnorm(1,0,50))
>> leaves <- replicate(500, releaf(leaf), simplify=FALSE)
>>
>> # Make into gpc.poly class:
>> leaves <- lapply(leaves, as, "gpc.poly")
>>
>> # Make union .
>> system.time({
>> leavesoutline <- union(leaves[[1]], leaves[[2]])
>> for(i in 3:length(leaves))leavesoutline <- union(leavesoutline,
>> leaves[[i]])
>> })
>> # about 1sec here.
>>
>> # Check it:
>> plot(leavesoutline)
>>
>>
>>
>> thanks!
>>
>> Remko
>>
>>
>> -
>> Remko Duursma
>> Research Lecturer
>>
>> Centre for Plants and the Environment
>> University of Western Sydney
>> Hawkesbury Campus
>> Richmond NSW 2753
>>
>> Dept of Biological Science
>> Macquarie University
>> North Ryde NSW 2109
>> Australia
>>
>> Mobile: +61 (0)422 096908
>> www.remkoduursma.com
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
>
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Faster union of polygons?
Reduce might work. Not sure about the speed advantages though. It does
simplify code.
Unionall <- function(x) Reduce('union', x)
leaveout <- Unionall(leaves)
On Tue, Jun 1, 2010 at 9:53 PM, Remko Duursma wrote:
> Dear R-helpers,
>
> thanks for yesterday's speeding-up tip. Here is my next query:
>
> I have lots of polygons (not necessarily convex ones, and they never
> have holes) given by x,y coordinates.
>
> I want to get the polygon that is the union of these polygons. This is
> my current method, but I am hoping there is a faster method (up to
> thousands of polygons, each with ca. 40 xy points).
>
> Example:
>
> library(gpclib)
>
> # A polygon
> leaf <- structure(c(0, 1, 12.9, 16.5, 18.8, 17, 16.8, 15.5, 12.1, 8.2,
> 6.3, 5, 2, 0, -1.5, -4.3, -6.6, -10.3, -14.8, -19.4, -22.2, -23.5,
> -22.2, -17.6, -7.8, 0, 0, -2.4, 2.8, 8.9, 19.9, 33.9, 34.8, 40.4,
> 49.7, 69.2, 77.4, 83.4, 91.4, 99, 92.8, 87.3, 81.2, 71.1, 57.6,
> 45.4, 39.2, 26, 15.6, 5.3, 0.6, 0), .Dim = c(26L, 2L), .Dimnames = list(
>NULL, c("X", "Y")))
>
> # Lots of polygons:
> releaf <-
> function(leaf)cbind(leaf[,1]+rnorm(1,0,50),leaf[,2]+rnorm(1,0,50))
> leaves <- replicate(500, releaf(leaf), simplify=FALSE)
>
> # Make into gpc.poly class:
> leaves <- lapply(leaves, as, "gpc.poly")
>
> # Make union .
> system.time({
> leavesoutline <- union(leaves[[1]], leaves[[2]])
> for(i in 3:length(leaves))leavesoutline <- union(leavesoutline,
> leaves[[i]])
> })
> # about 1sec here.
>
> # Check it:
> plot(leavesoutline)
>
>
>
> thanks!
>
> Remko
>
>
> -
> Remko Duursma
> Research Lecturer
>
> Centre for Plants and the Environment
> University of Western Sydney
> Hawkesbury Campus
> Richmond NSW 2753
>
> Dept of Biological Science
> Macquarie University
> North Ryde NSW 2109
> Australia
>
> Mobile: +61 (0)422 096908
> www.remkoduursma.com
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Faster union of polygons?
Baptiste,
thanks for the tip but this would give me an approximate union, and I
really need a (nearly) exact one. I am also not sure how to set the
alpha parameter in a non-arbitrary way.
Remko
-
Remko Duursma
Research Lecturer
Centre for Plants and the Environment
University of Western Sydney
Hawkesbury Campus
Richmond NSW 2753
Dept of Biological Science
Macquarie University
North Ryde NSW 2109
Australia
Mobile: +61 (0)422 096908
www.remkoduursma.com
On Wed, Jun 2, 2010 at 3:36 PM, baptiste auguie
wrote:
> Hi,
>
> I think you could use a concave hull from the alphahull package,
>
> http://yihui.name/en/2010/04/alphahull-an-r-package-for-alpha-convex-hull/
>
> It may be difficult to find the right parameters if the polygons
> differ widely in edge lengths, though.
>
>
> HTH,
>
> baptiste
>
> On 2 June 2010 03:53, Remko Duursma wrote:
>> Dear R-helpers,
>>
>> thanks for yesterday's speeding-up tip. Here is my next query:
>>
>> I have lots of polygons (not necessarily convex ones, and they never
>> have holes) given by x,y coordinates.
>>
>> I want to get the polygon that is the union of these polygons. This is
>> my current method, but I am hoping there is a faster method (up to
>> thousands of polygons, each with ca. 40 xy points).
>>
>> Example:
>>
>> library(gpclib)
>>
>> # A polygon
>> leaf <- structure(c(0, 1, 12.9, 16.5, 18.8, 17, 16.8, 15.5, 12.1, 8.2,
>> 6.3, 5, 2, 0, -1.5, -4.3, -6.6, -10.3, -14.8, -19.4, -22.2, -23.5,
>> -22.2, -17.6, -7.8, 0, 0, -2.4, 2.8, 8.9, 19.9, 33.9, 34.8, 40.4,
>> 49.7, 69.2, 77.4, 83.4, 91.4, 99, 92.8, 87.3, 81.2, 71.1, 57.6,
>> 45.4, 39.2, 26, 15.6, 5.3, 0.6, 0), .Dim = c(26L, 2L), .Dimnames = list(
>> NULL, c("X", "Y")))
>>
>> # Lots of polygons:
>> releaf <- function(leaf)cbind(leaf[,1]+rnorm(1,0,50),leaf[,2]+rnorm(1,0,50))
>> leaves <- replicate(500, releaf(leaf), simplify=FALSE)
>>
>> # Make into gpc.poly class:
>> leaves <- lapply(leaves, as, "gpc.poly")
>>
>> # Make union .
>> system.time({
>> leavesoutline <- union(leaves[[1]], leaves[[2]])
>> for(i in 3:length(leaves))leavesoutline <- union(leavesoutline, leaves[[i]])
>> })
>> # about 1sec here.
>>
>> # Check it:
>> plot(leavesoutline)
>>
>>
>>
>> thanks!
>>
>> Remko
>>
>>
>> -
>> Remko Duursma
>> Research Lecturer
>>
>> Centre for Plants and the Environment
>> University of Western Sydney
>> Hawkesbury Campus
>> Richmond NSW 2753
>>
>> Dept of Biological Science
>> Macquarie University
>> North Ryde NSW 2109
>> Australia
>>
>> Mobile: +61 (0)422 096908
>> www.remkoduursma.com
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Faster union of polygons?
Hi,
I think you could use a concave hull from the alphahull package,
http://yihui.name/en/2010/04/alphahull-an-r-package-for-alpha-convex-hull/
It may be difficult to find the right parameters if the polygons
differ widely in edge lengths, though.
HTH,
baptiste
On 2 June 2010 03:53, Remko Duursma wrote:
> Dear R-helpers,
>
> thanks for yesterday's speeding-up tip. Here is my next query:
>
> I have lots of polygons (not necessarily convex ones, and they never
> have holes) given by x,y coordinates.
>
> I want to get the polygon that is the union of these polygons. This is
> my current method, but I am hoping there is a faster method (up to
> thousands of polygons, each with ca. 40 xy points).
>
> Example:
>
> library(gpclib)
>
> # A polygon
> leaf <- structure(c(0, 1, 12.9, 16.5, 18.8, 17, 16.8, 15.5, 12.1, 8.2,
> 6.3, 5, 2, 0, -1.5, -4.3, -6.6, -10.3, -14.8, -19.4, -22.2, -23.5,
> -22.2, -17.6, -7.8, 0, 0, -2.4, 2.8, 8.9, 19.9, 33.9, 34.8, 40.4,
> 49.7, 69.2, 77.4, 83.4, 91.4, 99, 92.8, 87.3, 81.2, 71.1, 57.6,
> 45.4, 39.2, 26, 15.6, 5.3, 0.6, 0), .Dim = c(26L, 2L), .Dimnames = list(
> NULL, c("X", "Y")))
>
> # Lots of polygons:
> releaf <- function(leaf)cbind(leaf[,1]+rnorm(1,0,50),leaf[,2]+rnorm(1,0,50))
> leaves <- replicate(500, releaf(leaf), simplify=FALSE)
>
> # Make into gpc.poly class:
> leaves <- lapply(leaves, as, "gpc.poly")
>
> # Make union .
> system.time({
> leavesoutline <- union(leaves[[1]], leaves[[2]])
> for(i in 3:length(leaves))leavesoutline <- union(leavesoutline, leaves[[i]])
> })
> # about 1sec here.
>
> # Check it:
> plot(leavesoutline)
>
>
>
> thanks!
>
> Remko
>
>
> -
> Remko Duursma
> Research Lecturer
>
> Centre for Plants and the Environment
> University of Western Sydney
> Hawkesbury Campus
> Richmond NSW 2753
>
> Dept of Biological Science
> Macquarie University
> North Ryde NSW 2109
> Australia
>
> Mobile: +61 (0)422 096908
> www.remkoduursma.com
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
