Re: [R] Faster union of polygons?

From: Remko Duursma <remkoduursma_at_gmail.com>
Date: Fri, 04 Jun 2010 16:01:48 +1000

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 <mtmorgan_at_fhcrc.org> 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 <nikhil.list_at_gmail.com> 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 <remkoduursma_at_gmail.com>
>>> 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_at_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_at_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_at_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. Received on Fri 04 Jun 2010 - 06:05:12 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 04 Jun 2010 - 06:20:28 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive