sorry... my attached image had wrong annotations. The correctly annotated image is attached to this email.
On Fri, Sep 10, 2010 at 11:27 AM, P Kishor <[email protected]> wrote: > On Fri, Sep 10, 2010 at 2:44 AM, Matthew Kenworthy > <[email protected]> wrote: >> >>> I wrote a PDL vectorized version of the 'points in a polygon' program a >>> few years ago for my own use. [BTW Matt the problem was selecting points in >>> a area of a color-color diagram] >>> >>> Here is the code. i/o are all vectors. Read it and weep :-) >>> >>> Karl >>> >>> # Test points tx,ty to be in a polygon following >>> # http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html >> >> *jaw bounces off floor* Okay, that's one slick algorithm. Puneet, does this >> do the job you're after? >> > > > Well, no need to be so modest, Matt, your solution stacks up quite > well against Karl's. So, first, the problem -- > > Given a 2D piddle and a set of coordinates that make an irregular poly > in that piddle, return a mask piddle with all elements within the > irregular poly set to 1, and all outside set to 0. > > I tested four programs -- Matt's rather ingenious solution, Karl's > vectorized point-in-poly, my own non-PDL naive point-in-poly, and my > not-so-naive non-PDL point-in-poly. > > The naive version calculates a point-in-poly for every point in the 2D > input piddle. > > The not-so-naive version first calculates a bounding box for the > irregular poly, and then calculates point-in-poly for every point in > the bounding box. I lifted point-in-poly and bounding-box algorithms > from the "Mastering Algorithms with Perl" book. > > I benchmarked the four programs with a 1000 x 1000 input piddle and > small convex and concave polygons (really small, only 5-6 points > each). > > First, the results -- > > In terms of accuracy, although no solution is absolutely accurate (per > my definition of accuracy), Matt's solution seems to provide better > accuracy. The point-in-poly approach seems to get tripped by and > ignore vertical and horizontal lines. Karl's and my both naive and > not-so solutions produce identical results as they are all > point-in-poly approaches. > > See the attached images. The black lines make the irregular poly, the > dark orange are the selected elements, the light orange elements are > those that should have also been selected but were skipped. > > That said, Matt's solution, admittedly, works only for convex hulls. > The pnp solutions work on concave hulls as well. > > In terms of speed, Matt's and Karl's solutions, and, surprise (!), my > not-so-naive program, are all almost equally fast. Actually, Matt's > solution is twice as fast as Karl's. I ran Devel::NYTProf, and > according to that, Matt's solution takes 348 ms, while Karl's takes > 962ms. > > 'karl' took : 1 wallclock secs ( 0.49 usr + 0.57 sys = 1.06 CPU) > 'matt' took : 0 wallclock secs ( 0.25 usr + 0.05 sys = 0.30 CPU) > 'naive' took : 16 wallclock secs (14.36 usr + 0.06 sys = 14.42 CPU) > 'not_so' took : 0 wallclock secs ( 0.00 usr + 0.00 sys = 0.00 CPU) > > > Most of the time taken by Matt's solution is spent in creating the > mask piddle, 150 ms out of the total 348 ms. The culprit is > > my $mask = ( > (yvals($pdl)->dummy(0)) > ($m * (xvals($pdl)->dummy(0)) + $c) > ); > > Most of the time taken by Karl's solution is spent in the sumover > routine, 802 ms out of the total 962 ms. The offending line is > > my $c = sumover( > (($verty > $testy) != ($vertyj > $testy)) & > ($testx < ($vertxj - $vertx) * ($testy - $verty) / > ($vertyj - $verty) + $vertx) > ) % 2; > > My naive, non-PDL point-in-poly is the slowest, as expected. Making > 1000_000 loops is just silly. > > However, my not_so_naive, non-PDL point-in-poly is almost as fast as > the PDL versions! Restricting the loops to the bounding box of the > irregular poly reduces the loops to 30, making non-PDL solution to be > just as fast as PDL. Of course, in reality, the irregular poly is > going to be much bigger. But still, even at a bounding box of 100 x > 100, the non-PDL solution is very fast, under a second. Might be worth > it for someone who doesn't have access to PDL, or doesn't have > Matt/Karl to help out. > > Now, what would be really neat is if an R-Tree implementation was > added to PDL. Someone surprise me and tell me it has already been > done. > > > Kudos to Matt for teaching me step by step how his algorithm works... > it is really very smart (or, I am easily impressed). Kudos to Karl, as > usual, for coming out of nowhere with a jaw-dropper. > > > > > > > -- > Puneet Kishor http://www.punkish.org > Carbon Model http://carbonmodel.org > Charter Member, Open Source Geospatial Foundation http://www.osgeo.org > Science Commons Fellow, http://sciencecommons.org/about/whoweare/kishor > Nelson Institute, UW-Madison http://www.nelson.wisc.edu > ----------------------------------------------------------------------- > Assertions are politics; backing up assertions with evidence is science > ======================================================================= > -- Puneet Kishor http://www.punkish.org Carbon Model http://carbonmodel.org Charter Member, Open Source Geospatial Foundation http://www.osgeo.org Science Commons Fellow, http://sciencecommons.org/about/whoweare/kishor Nelson Institute, UW-Madison http://www.nelson.wisc.edu ----------------------------------------------------------------------- Assertions are politics; backing up assertions with evidence is science =======================================================================
<<attachment: irr_poly_select.png>>
_______________________________________________ Perldl mailing list [email protected] http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
