Re: [R] Faster union of polygons?

From: Martin Morgan <mtmorgan_at_fhcrc.org>
Date: Thu, 03 Jun 2010 19:17:34 -0700

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 - 02:17:57 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:10: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