Hi, Chris.
I updated the documentation for polyfill and polyfillv, added the Perl
subroutine PDL::polyfill, and added the pnpoly option to polyfill and
polyfillv. I have attached the git diff file generated against the latest
repository. It builds and tests without issue against Perl 5.14.2 on Ubuntu
12.04 x86_64.
Many thanks.
- Tim
On Thu, Jun 14, 2012 at 2:50 PM, Chris Marshall <[email protected]>wrote:
> Thanks for the performance measurements to
> confirm our hypothesis.
>
> As for variability in the results, the key is
> that pnpoly checks points against an ideal
> polygon described by vertex coordinates
> and *not* pixels against the rendered image
> of the polygon.
>
> I think updating the documentation is sufficient.
> We could also modify the polyfill routine to have
> a choice of engines based on an option setting
> or flag.
>
> --Chris
>
> On Thu, Jun 14, 2012 at 1:30 PM, Tim Haines <[email protected]>
> wrote:
> > Hi, Chris.
> >
> >>> I suggest changing the docs for polyfill and polyfillv to note that the
> >>> algorithm is
> >>> faster/less memory intensive (is it?) than pnpoly but slightly less
> >>> accurate.
> >
> > Using Benchmark::cmpthese, I found the following with 1e4 iterations on
> the
> > triangle example:
> >
> > Rate pnpoly polyfill
> > pnpoly 5263/s -- -95%
> > polyfill 111111/s 2011% --
> >
> > Indeed, pnpoly is quite a bit slower than polyfill and uses more memory
> due
> > to the intermediates xvals/yvals and the sumover() call made internally
> by
> > pnpoly.
> >
> >>> I ran your test code and it died after the following output...
> >
> > Are you using Perl 5.16? I forgot to mention that I am using Perl 5.14.2
> and
> > PDL 2.4.10.
> >
> >
> >>> It appears that polyfillv is not marked as an lvalue sub.
> >
> > Thanks for updating that.
> >
> >
> >>> The differences between the two results appear to be whether or not the
> >>> edge pixels are included and the pnpoly does
> >>> appear to be "fully" correct in generating the set of interior points
> to
> >>> the polygon.
> >
> > For the most part, yes. I have attached a marked-up version of the output
> > highlighting the different vertices and edges marked in red and green,
> > respectively. However, it seems the problem is even worse than I had
> > originally discovered. After a more careful inspection, I found that
> > polyfill gives the correct answer for the area of the triangle, but the
> > reverse is true for the square!
> >
> > Testing the triangle (actual area = 0.5*12*6 = 36)
> > pnpoly area = 42
> > polyfill area = 36
> >
> > Testing the square (actual area = 4*4=16)
> > pnpoly area = 16
> > polyfill area = 20
> >
> > Given that neither method is consistently reliable on these simple
> > geometries, I tested the original pnpoly implementation that Karl
> vectorized
> > to create PDL::pnpoly
> > (http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
> ).
> > Karl's version faithfully represents the original implementation. Hence,
> it
> > is an issue inherent to the method being employed. My guess would be it
> is
> > the inclusion of the "Simulation of Simplicity" that is causing the
> > preferential selection of the horizontal edge in the triangle. To test
> this,
> > I performed the test again with a rotated triangle that has no
> horizontal or
> > vertical edges ($px = pdl(3,9,13) and $py = pdl(3,10,5)).
> >
> > Test the rotated triangle (actual area = ~34)
> > pnpoly area = 28
> > polyfill area = 37
> >
> > This result is also congruent with the original pnpoly implementation,
> and,
> > unsurprisingly, the "discretized" area isn't quite right. It seems that
> the
> > "Simulation of Simplicity" is only looking "into" the polygon to
> determine
> > if it is on a boundary. This would explain the inclusion of the
> horizontal
> > edge of the regular triangle and the "top" of the square, but not
> including
> > the "bottom" of the square.
> >
> > I am uncertain how to proceed. polyfill is consistent, but incorrect.
> pnpoly
> > is inconsistent, but mostly correct. Neither of these is a very good
> > position to be in. However, I am tending to agree with your original
> > assessment of changing the documentation to state that polyfill[v] does
> > include *some* of the vertices and edges.
> >
> > Many thanks.
> >
> > - Tim
> >
> >
> >
> >
> >
> > On Thu, Jun 14, 2012 at 7:22 AM, Chris Marshall <[email protected]>
> > wrote:
> >>
> >> I ran your test code and it died after
> >> the following output:
> >>
> >> Testing the triangle
> >>
> >> [
> >> [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
> >> [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
> >> [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
> >> [0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0]
> >> [0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0]
> >> [0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0]
> >> [0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0]
> >> [0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0]
> >> [0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0]
> >> [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
> >> [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
> >> ]
> >> [0 0 0 12 10 8 6 4 2 0 0]
> >> Can't modify non-lvalue subroutine call at polyfill-test.pl line 31.
> >>
> >> It appears that polyfillv is not marked as an
> >> lvalue sub. If not, I would consider that a
> >> bug. The differences between the two
> >> results appear to be whether or not the edge
> >> pixels are included and the pnpoly does
> >> appear to be "fully" correct in generating the
> >> set of interior points to the polygon.
> >>
> >> --Chris
> >>
> >> On Thu, Jun 14, 2012 at 8:08 AM, Chris Marshall <[email protected]
> >
> >> wrote:
> >> > Putting polyfill and polyfillv on the path to
> >> > deprecate the old implementation and move
> >> > to the pnpoly one could be done as well.
> >> >
> >> > --Chris
> >> >
> >> > On Wed, Jun 13, 2012 at 5:50 PM, Chris Marshall <
> [email protected]>
> >> > wrote:
> >> >> Hi Tim-
> >> >>
> >> >> I note that polyfill and polyfillv have been part of PDL
> >> >> since 2000 while pnpoly was only added Dec 2011.
> >> >>
> >> >> I'm a bit leery of changing something that has been
> >> >> part of the PDL distribution for such a long time without
> >> >> a clear need to do so. I suggest changing the docs
> >> >> for polyfill and polyfillv to note that the algorithm is
> >> >> faster/less memory intensive (is it?) than pnpoly but
> >> >> slightly less accurate.
> >> >>
> >> >> Any other opinions on this? BTW, I assume that the
> >> >> difference in on the edges of the regions and no a
> >> >> big difference between the pnpoly result.
> >> >>
> >> >> --Chris
> >> >>
> >> >> On Wed, Jun 13, 2012 at 3:45 PM, Tim Haines
> >> >> <[email protected]> wrote:
> >> >>> Greetings, all.
> >> >>>
> >> >>> I was playing around with Image2D::pnpoly and Image2D::polyfill (and
> >> >>> polyfillv), and I noticed that the latter two do not have the same
> >> >>> behavior
> >> >>> as pnpoly. I have some tests below.
> >> >>>
> >> >>> use strict;
> >> >>> use warnings;
> >> >>> use PDL;
> >> >>> use PDL::Image2D;
> >> >>>
> >> >>> # A triangle
> >> >>> my $px = pdl(3,15,9);
> >> >>> my $py = pdl(3,3,9);
> >> >>> print "Testing the triangle\n";
> >> >>> &testPoly($px,$py);
> >> >>>
> >> >>> # A square (be careful to draw it convexly!)
> >> >>> $px = pdl(4,4,8,8);
> >> >>> $py = pdl(4,8,8,4);
> >> >>> print "Testing the square\n";
> >> >>> &testPoly($px,$py);
> >> >>>
> >> >>> sub testPoly
> >> >>> {
> >> >>> my $px = shift;
> >> >>> my $py = shift;
> >> >>> my $points = $px->cat($py)->xchg(0,1);
> >> >>>
> >> >>> my $img1 = zeroes(20, 11);
> >> >>> my $img2 = zeroes(20, 11);
> >> >>> my $img3 = zeroes(20, 11);
> >> >>> my $img4 = zeroes(20, 11);
> >> >>>
> >> >>> my $img5 = pnpoly($img1->xvals,$img1->yvals,$px,$py);
> >> >>> print $img5, $img5->sumover, "\n";
> >> >>> $img2->polyfillv($points) .= 1;
> >> >>> print $img2, $img2->sumover, "\n";
> >> >>> $img3->polyfill($points, 1);
> >> >>> print $img3, $img3->sumover, "\n";
> >> >>> }
> >> >>>
> >> >>>
> >> >>> From the output, I believe that pnpoly produces the correct answer
> >> >>> insofar
> >> >>> as it produces the correct area for the square. Moreover, the
> >> >>> documentation
> >> >>> for polyfill and polyfillv states that "fill the area inside the
> given
> >> >>> polygon with a given colour" and "return the (dataflown) area of an
> >> >>> image
> >> >>> within a polygon," respectively. However, the area selected for
> color
> >> >>> filling in polyfill contains more than the points interior to the
> >> >>> polygon,
> >> >>> and similarly for polyfillv.
> >> >>>
> >> >>> After some investigation into PDL/Lib/Image2D, I found that
> polyfillv
> >> >>> uses
> >> >>> polyfill, so their identical behavior is expected.
> >> >>>
> >> >>> My proposed fix is to use pnpoly to determine the interior points.
> >> >>> Then
> >> >>> polyfill becomes simply
> >> >>>
> >> >>> sub PDL::polyfill
> >> >>> {
> >> >>> my ($im, $ps, $col) = @_;
> >> >>> barf("ps must contain pairwise points.\n") unless $ps->getdim(0)
> >> >>> == 2;
> >> >>> barf("color not specified.\n") unless $col;
> >> >>> $im->where($im->pnpoly($im->xvals, $im->yvals,
> $ps->slice('0,:'),
> >> >>> $ps->slice('1,:'))) .= $col;
> >> >>> }
> >> >>>
> >> >>> and polyfillv becomes
> >> >>>
> >> >>> sub PDL::polyfillv
> >> >>> {
> >> >>> my ($im, $ps) = @_;
> >> >>> barf("ps must contain pairwise points.\n") unless $ps->getdim(0)
> >> >>> == 2;
> >> >>> return $im->where($im->pnpoly($im->xvals, $im->yvals,
> >> >>> $ps->slice('0,:'),
> >> >>> $ps->slice('1,:')));
> >> >>> }
> >> >>>
> >> >>> This will remove the need for the PP version of polyfill. Assuredly,
> >> >>> these
> >> >>> fixes are not optimal. If $im is large, then two large intermediates
> >> >>> ($im->xvals and $im->yvals) are created.
> >> >>>
> >> >>> I did not file this as a bug report, I wanted to run it by everyone
> to
> >> >>> see
> >> >>> if I was missing something. Additionally, I am not closely familiar
> >> >>> with the
> >> >>> coding standards, so I would appreciate any tips on that, as well.
> >> >>>
> >> >>> Many thanks.
> >> >>>
> >> >>> - Tim
> >> >>>
> >> >>>
> >> >>>
> >> >>>
> >> >>> _______________________________________________
> >> >>> Perldl mailing list
> >> >>> [email protected]
> >> >>> http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
> >> >>>
> >
> >
>
diff --git a/Lib/Image2D/image2d.pd b/Lib/Image2D/image2d.pd
index 3f6318e..f6269f3 100644
--- a/Lib/Image2D/image2d.pd
+++ b/Lib/Image2D/image2d.pd
@@ -1193,7 +1193,7 @@ pp_addpm(<<'EOPM');
$mask = pnpoly($x, $y, $px, $py);
-For a closed polygon determined by the sequence of points in {$py,$py}
+For a closed polygon determined by the sequence of points in {$px,$py}
the output of pnpoly is a mask corresponding to whether or not each
coordinate (x,y) in the set of test points, {$x,$y}, is in the interior
of the polygon. This is the 'points in a polygon' algorithm from
@@ -1258,12 +1258,41 @@ pp_def('polyfill',
=for ref
-fill the area inside the given polygon with a given colour
+fill the area of the given polygon with a given colour.
This function works inplace, i.e. modifies C<im>.
+=for usage
+
+ polyfill($im,$ps,$colour,[pnpoly]);
+
+The method of determining which points lie inside of the polygon used here
+is not as strict as the method used in L<pnpoly|pnpoly>. Often, it includes
vertices
+and edge points. However, this method is faster, has a smaller memory
+footprint, and works inplace. If you wish to use the C<pnpoly> method,
+set pnpoly => 1.
+
+=for example
+
+ # Make a convex 3x3 square of 1s in an image using the pnpoly method
+ $ps = pdl([3,3],[3,6],[6,6],[6,3]);
+ polyfill($im,$ps,1,'pnpoly'=>1);
+
=cut
+sub PDL::polyfill {
+ my ($im,$ps,$val,%pnpoly) = @_;
+ if($pnpoly{'pnpoly'}) {
+ my $msk =
pnpoly($im->xvals,$im->yvals,$ps->slice('0,:'),$ps->slice('1,:'));
+ $im->where($msk) .= $val;
+ } else {
+ polyfill($im,$ps,$val);
+ }
+ return $im;
+}
+
+*polyfill = \&PDL::polyfill;
+
EOD
);
@@ -1274,21 +1303,27 @@ pp_addpm(<<'EOPM');
=for ref
-return the (dataflown) area of an image within a polygon
+return the (dataflown) area of an image described by a polygon
+
+=for usage
+
+ $im->polyfillv($ps,[pnpoly]);
+
+See C<polyfill> for details on usage.
=for example
# increment intensity in area bounded by $poly
$im->polyfillv($pol)++; # legal in perl >= 5.6
- # compute average intensity within area bounded by $poly
+ # compute average intensity of area bounded by $poly
$av = $im->polyfillv($poly)->avg;
=cut
sub PDL::polyfillv :lvalue {
- my ($im, $ps) = @_;
+ my ($im, $ps, %pnpoly) = @_;
my $msk = zeroes(long,$im->dims);
- polyfill($msk, $ps, 1);
+ PDL::polyfill($msk, $ps, 1, %pnpoly);
return $im->where($msk == 1);
}
*polyfillv = \&PDL::polyfillv;
_______________________________________________
Perldl mailing list
[email protected]
http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
