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
=======================================================================
<<attachment: irr_poly_select.png>>
_______________________________________________ Perldl mailing list [email protected] http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
