Skip to content

Commit

Permalink
merge in pnpoly algorithm into PDL::Image2D
Browse files Browse the repository at this point in the history 
This commit represents the work by Tim Haines to
implement a PP version of the pnpoly algorithm and
to add that as an optional engine for the polyfill
and polyfillv routines.
  • Loading branch information
@devel-chm
Tim Haines authored and devel-chm committed Jul 22, 2012
1 parent 3cff6f1 commit 335023d
Show file tree
Hide file tree
Showing 3 changed files with 229 additions and 45 deletions.
24 changes: 24 additions & 0 deletions Basic/Pod/PP.pod
View file
Expand Up @@ -1635,6 +1635,30 @@ two. To use this, simply have the following line at some point in your
However, don't use this if you use L<Module::Build::PDL>. See that
module's documentation for details.

=head1 Making your PP function "private"

Let's say that you have a function in your module called PDL::foo that
uses the PP function C<bar_pp> to do the heavy lifting. But you don't
want to advertise that C<bar_pp> exists. To do this, you must move your
PP function to the top of your module file, then call

pp_export_nothing()

to clear the C<EXPORT> list. To ensure that no documentation (even the
default PP docs) is generated, set

Doc => undef

and to prevent the function from being added to the symbol table, set

PMFunc => ''

in your pp_def declaration (see Image2D.pd for an example). This will
effectively make your PP function "private." However, it is I<always>
accessible via PDL::bar_pp due to Perl's module design. But making
it private will cause the user to go very far out of his or her way
to use it, so he or she shoulders the consequences!

=head1 Slice operation

The slice operation section of this manual is provided using
Expand Down
230 changes: 186 additions & 44 deletions Lib/Image2D/image2d.pd
View file
Expand Up @@ -44,6 +44,67 @@ the copyright notice should be included in the file.

EOD

#################################################
# BEGIN INTERNAL FUNCTION DECLARATIONS #
#################################################
pp_def('polyfill_pp',
HandleBad => 0, # a marker
Pars => 'int [o,nc] im(m,n); float ps(two=2,np); int col()',
Code => 'int ierr = 0, nerr;
threadloop %{
polyfill($P(im), $SIZE(m), $SIZE(n), $P(ps), $SIZE(np), $col(), &nerr);
ierr = ierr < nerr ? nerr : ierr;
%}
if (ierr) warn("errors during polygonfilling");
',
Doc => undef,
PMFunc => ''
);

my %pnpolyFields = (
'pnpoly_pp' => {'pars' => 'a(m,n); ps(k,l); int [o] msk(m,n)', 'special' => '$msk() = c;'},
'pnpolyfill_pp' => {'pars' => '[o,nc] a(m,n); ps(k,l); int col()', 'special' => 'if(c) { $a() = $col(); }'}
);

for my $name (keys %pnpolyFields) {
pp_def($name,
HandleBad => 0,
PMFunc => '',
Doc => undef,
Pars => $pnpolyFields{$name}->{'pars'},
Code => '
int i, j, c, nvert;
nvert = $SIZE(l);

#define VERTX(q) $ps(k=>0,l=>q)
#define VERTY(q) $ps(k=>1,l=>q)

threadloop %{
loop(n) %{
loop(m) %{
c = 0;
for(i=0,j=nvert-1;i<nvert;j=i++) {
if( ((VERTY(i)>n) != (VERTY(j)>n)) &&
(m < (VERTX(j)-VERTX(i)) * (n-VERTY(i)) / (VERTY(j)-VERTY(i)) + VERTX(i)) )
c = !c;
}
' . $pnpolyFields{$name}->{'special'} .'
%}
%}
%}

#undef VERTX
#undef VERTY
'
);
}

pp_export_nothing(); # Clear the export list

#################################################
# END INTERNAL FUNCTION DECLARATIONS #
#################################################


pp_addhdr('

Expand Down Expand Up @@ -1180,6 +1241,63 @@ void polyfill(PDL_Long *image, int wx, int wy, float *ps, int n,

');

pp_add_exported('polyfill');
pp_addpm(<<'EOPM');
=head1 polyfill

=for ref

fill the area of the given polygon with the given colour.

This function works inplace, i.e. modifies C<im>.

=for usage

polyfill($im,$ps,$colour,[\%options]);

The default method of determining which points lie inside of the polygon used
is not as strict as the method used in L<pnpoly|pnpoly>. Often, it includes vertices
and edge points. Set the C<Method> option to change this behaviour.

=for option

Method - Set the method used to determine which points lie in the polygon.
=> Default - internal PDL algorithm
=> pnpoly - use the L<pnpoly|pnpoly> algorithm

=for example

# Make a convex 3x3 square of 1s in an image using the pnpoly algorithm
$ps = pdl([3,3],[3,6],[6,6],[6,3]);
polyfill($im,$ps,1,{'Method' =>'pnpoly'});

=cut
sub PDL::polyfill {
my $opt;
$opt = pop @_ if ref($_[-1]) eq 'HASH';
barf('Usage: polyfill($im,$ps,$colour,[\%options])') unless @_==3;
my ($im,$ps,$colour) = @_;

if($opt) {
use PDL::Options qw();
my $parsed = PDL::Options->new({'Method' => undef});
$parsed->options($opt);
if( $parsed->current->{'Method'} eq 'pnpoly' ) {
PDL::pnpolyfill_pp($im,$ps,$colour);
}
}
else
{
PDL::polyfill_pp($im,$ps,$colour);
}
return $im;

}

*polyfill = \&PDL::polyfill;

EOPM

pp_add_exported('', 'pnpoly');
pp_addpm(<<'EOPM');

Expand All @@ -1190,10 +1308,13 @@ pp_addpm(<<'EOPM');
'points in a polygon' selection from a 2-D piddle

=for usage


$mask = $img->pnpoly($ps);

# Old style, do not use
$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
Expand All @@ -1203,10 +1324,16 @@ and vectorized for PDL by Karl Glazebrook.
=for example

# define a 3-sided polygon (a triangle)
$ps = pdl([3, 3], [20, 20], [34, 3]);

# $tri is 0 everywhere except for points in polygon interior
$tri = $img->pnpoly($ps);

With the second form, the x and y coordinates must also be specified.
B< I<THIS IS MAINTAINED FOR BACKWARD COMPATIBILITY ONLY> >.

$px = pdl( 3, 20, 34 );
$py = pdl( 3, 20, 3 );

$img = zeros(40, 40); # create test image
$x = $img->xvals; # get x pixel coords
$y = $img->yvals; # get y pixel coords

Expand All @@ -1229,67 +1356,82 @@ and vectorized for PDL by Karl Glazebrook.
# operator does not "short circuit".

sub PDL::pnpoly {
my ($tx, $ty, $vertx, $verty) = @_;
my $testx = PDL::Core::topdl($tx)->dummy(0);
my $testy = PDL::Core::topdl($ty)->dummy(0);
my $vertxj = PDL::Core::topdl($vertx)->rotate(1);
my $vertyj = PDL::Core::topdl($verty)->rotate(1);
my $mask = ( ($verty>$testy) == ($vertyj>$testy) );
my $c = sumover( ! $mask & ( $testx < ($vertxj-$vertx) * ($testy-$verty)
/ ($vertyj-$verty+$mask) + $vertx) ) % 2;
return $c;
barf('Usage: $mask = pnpoly($img, $ps);') unless(@_==2 || @_==4);
my ($tx, $ty, $vertx, $verty) = @_;

# if only two inputs, use the pure PP version
unless( defined $vertx ) {
barf("ps must contain pairwise points.\n") unless $ty->getdim(0) == 2;

# Input mapping: $img => $tx, $ps => $ty
return PDL::pnpoly_pp($tx,$ty);
}

my $testx = PDL::Core::topdl($tx)->dummy(0);
my $testy = PDL::Core::topdl($ty)->dummy(0);
my $vertxj = PDL::Core::topdl($vertx)->rotate(1);
my $vertyj = PDL::Core::topdl($verty)->rotate(1);
my $mask = ( ($verty>$testy) == ($vertyj>$testy) );
my $c = sumover( ! $mask & ( $testx < ($vertxj-$vertx) * ($testy-$verty)
/ ($vertyj-$verty+$mask) + $vertx) ) % 2;
return $c;
}

*pnpoly = \&PDL::pnpoly;

EOPM

pp_def('polyfill',
HandleBad => 0, # a marker
Pars => 'int [o,nc] im(m,n); float ps(two=2,np); int col()',
Code => 'int ierr = 0, nerr;
threadloop %{
polyfill($P(im), $SIZE(m), $SIZE(n), $P(ps), $SIZE(np), $col(), &nerr);
ierr = ierr < nerr ? nerr : ierr;
%}
if (ierr) warn("errors during polygonfilling");
',
Doc => << 'EOD',

=for ref

fill the area inside the given polygon with a given colour

This function works inplace, i.e. modifies C<im>.

=cut

EOD
);

pp_add_exported('', 'polyfillv');
pp_addpm(<<'EOPM');

=head2 polyfillv

=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

polyfillv($im,$ps,[\%options]);

The default method of determining which points lie inside of the polygon used
is not as strict as the method used in L<pnpoly|pnpoly>. Often, it includes vertices
and edge points. Set the C<Method> option to change this behaviour.

=for option

Method - Set the method used to determine which points lie in the polygon.
=> Default - internal PDL algorithm
=> pnpoly - use the L<pnpoly|pnpoly> algorithm

=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
# increment intensity in area bounded by $poly using the pnpoly algorithm
$im->polyfillv($poly,{'Method'=>'pnpoly'})++; # legal in perl >= 5.6

# compute average intensity within area bounded by $poly using the default algorithm
$av = $im->polyfillv($poly)->avg;

=cut

sub PDL::polyfillv :lvalue {
my ($im, $ps) = @_;
my $msk = zeroes(long,$im->dims);
polyfill($msk, $ps, 1);
return $im->where($msk == 1);
my $opt;
$opt = pop @_ if ref($_[-1]) eq 'HASH';
barf('Usage: polyfillv($im,$ps,[\%options])') unless @_==2;

my ($im,$ps) = @_;
barf("ps must contain pairwise points.\n") unless $ps->getdim(0) == 2;

if($opt) {
use PDL::Options qw();
my $parsed = PDL::Options->new({'Method' => undef});
$parsed->options($opt);
return $im->where(PDL::pnpoly_pp($im, $ps)) if $parsed->current->{'Method'} eq 'pnpoly';
}

my $msk = zeroes(long,$im->dims);
PDL::polyfill_pp($msk, $ps, 1);
return $im->where($msk);
}
*polyfillv = \&PDL::polyfillv;

Expand Down
20 changes: 19 additions & 1 deletion t/image2d.t
View file
Expand Up @@ -3,7 +3,7 @@

use Test;
BEGIN {
plan tests => 22;
plan tests => 25;
}

use PDL;
Expand Down Expand Up @@ -147,9 +147,27 @@ ok($@ eq '');
my $px = pdl(0,3,1);
my $py = pdl(0,1,4);
my $im = zeros(5,5);
my $im2 = zeroes(5,5);
my $x = $im->xvals;
my $y = $im->yvals;
my $ps = $px->cat($py)->xchg(0,1);
my $im_mask = pnpoly($x,$y,$px,$py);
ok(sum($im_mask) == 5);
my $inpixels = pdl q[ 1 1 ; 1 2 ; 1 3 ; 2 1 ; 2 2 ];
ok(sum($inpixels - qsortvec(scalar whichND($im_mask))) == 0);

# Make sure the PDL pnpoly and the PP pnpoly give the same result
ok(all($im_mask == $im->pnpoly($ps)));

# Trivial test to make sure the polyfills using the pnpoly algorithm are working
$im .= 0;
polyfillv($im2,$ps,{'Method'=>'pnpoly'}) .= 22;
ok(all(polyfill($im,$ps,22,{'Method'=>'pnpoly'}) == $im2));


# Trivial test to make sure the polyfills are working
$im .= 0;
$im2 .= 0;
polyfillv($im2,$ps) .= 25;
polyfill($im,$ps,25);
ok(all($im == $im2));

0 comments on commit 335023d

Please sign in to comment.