File Coverage

lib/PDL/Image2D.pd

Criterion	Covered	Total	%
statement	8	8	100.0
branch	92	1104	8.3
condition			n/a
subroutine			n/a
pod			n/a
total	100	1112	8.9

use strict;
use warnings;
use PDL::Types qw(ppdefs_all);
my $A = [ppdefs_all];

{ no warnings 'once'; # pass info back to Makefile.PL
$PDL::Core::Dev::EXTRAS{$::PDLMOD}{OBJECT} .= join '', map " $::PDLBASE/$_\$(OBJ_EXT)", qw(rotate resample select misc);
$PDL::Core::Dev::EXTRAS{$::PDLMOD}{INC} .= qq{ "-I$::PDLBASE"};
}

pp_addpm({At=>'Top'},<<'EOD');
use strict;
use warnings;

=head1 NAME

PDL::Image2D - Miscellaneous 2D image processing functions

=head1 DESCRIPTION

Miscellaneous 2D image processing functions - for want
of anywhere else to put them.

=head1 SYNOPSIS

 use PDL::Image2D;

=cut

use PDL;  # ensure qsort routine available
use PDL::Math;
use Carp;

my %boundary2value = (Reflect=>1, Truncate=>2, Replicate=>3);
EOD

pp_addpm({At=>'Bot'},<<'EOD');

=head1 AUTHORS

Copyright (C) Karl Glazebrook 1997 with additions by Robin Williams
([e-mail]), Tim Jenness ([e-mail]),
and Doug Burke ([e-mail]).

All rights reserved. There is no warranty. You are allowed
to redistribute this software / documentation under certain
conditions. For details, see the file COPYING in the PDL
distribution. If this file is separated from the PDL distribution,
the copyright notice should be included in the file.

=cut

EOD

#################################################
#     BEGIN INTERNAL FUNCTION DECLARATIONS      #
#################################################

pp_def('polyfill_pp',
    HandleBad => 0, # a marker
    Pars => '[io] im(m,n); float ps(two=2,np); int col()',
    GenericTypes => ['L'],
    CHeader => qq{void polyfill(PDL_Long *image, int wx, int wy, float *ps, int n, PDL_Long col, int *ierr);\n},
    Code => 'int ierr = 0, nerr;
             broadcastloop %{
               polyfill($P(im), $SIZE(m), $SIZE(n), $P(ps), $SIZE(np), $col(), &nerr);
               ierr = PDLMAX(ierr, nerr);
             %}
             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' => '[io] a(m,n); ps(k,l); int col()', 'special' => 'if(c) { $a() = $col(); }'}
);

for my $name (sort keys %pnpolyFields) {
	pp_def($name,
		HandleBad => 0,
		PMFunc => '',
		Doc => undef,
		Pars => $pnpolyFields{$name}->{'pars'},
		Code => '
#define VERTX(q) $ps(k=>0,l=>q)
#define VERTY(q) $ps(k=>1,l=>q)

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

#undef VERTX
#undef VERTY
'
	);
}

pp_export_nothing();   # Clear the export list

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


pp_addhdr('

/* Fast Modulus with proper negative behaviour */

#define REALMOD(a,b) {while ((a)>=(b)) (a) -= (b); while ((a)<0) (a) += (b);}

#define X(symbol, ctype, ppsym, ...) \
  ctype quick_select_ ## ppsym(ctype arr[], int n);
PDL_TYPELIST_REAL(X)
#undef X
');

my %init =
    (
      i => { size => 'm_size', off => 'poff', init => '1-p_size' },
      j => { size => 'n_size', off => 'qoff', init => '1-q_size' },
      );

# requires 'int $var, ${var}2' to have been declared in the c code
# (along with [pq]off and [pq]_size)
#
sub init_map {
    my $var = shift;

    my $loop = $var;
    my $loop2 = "${var}2";

    my $href = $init{$var} ||
	die "ERROR: unknown variable sent to init_map()\n";
    my $size = $href->{size} ||
	die "ERROR: unable to find size for $var\n";
    my $off  = $href->{off} ||
	die "ERROR: unable to find off for $var\n";
    my $init = $href->{init} ||
	die "ERROR: unable to find init for $var\n";

    return
"for ( $loop = $init; $loop< $size; ${loop}++) {
    $loop2 = $loop + $off;
    switch (opt) {
       case 1:      /* REFLECT */
          if (${loop2}<0)
             $loop2 = -${loop2}-1;
          else if ($loop2 >= $size)
             $loop2 = 2*${size}-(${loop2}+1);
          break;
       case 2:      /* TRUNCATE */
          if (${loop2}<0 || ${loop2} >= $size)
             $loop2 = -1;
          break;
       case 3:      /* REPLICATE */
          if (${loop2}<0)
             $loop2 = 0;
          if (${loop2} >= $size)
             $loop2 = $size-1;
          break;
       default:
           REALMOD($loop2,$size);
    }
    map${var}\[$loop] = $loop2;
 }\n";

} # sub: init_map()

sub init_vars {
    my $sizevars = shift || [qw(m n p q)];
    my $compvars = shift || [];
    my $str = "int i,j, i2,j2, poff, qoff;";
    $str .=
	'int opt = $COMP(opt);
	' . join("\n", map "int ${_}_size = \$SIZE(${_});", @$sizevars) . '
	' . join("\n", map "int ${_}_size = \$COMP(${_}_size);", @$compvars) . '
         PDL_Indx *mapi = $P(mapi), *mapj = $P(mapj);
         if ((mapi==NULL) || (mapj==NULL)) $CROAK("Out of Memory");

         poff = p_size/2; mapi += p_size-1;
         qoff = q_size/2; mapj += q_size-1;
';

    return $str;
} # sub: init_vars()

pp_def('conv2d', Doc=><<'EOD',
=for ref

2D convolution of an array with a kernel (smoothing)

For large kernels, using a FFT routine,
such as L,
will be quicker.

=for usage

 $new = conv2d $old, $kernel, {OPTIONS}

=for example

 $smoothed = conv2d $image, ones(3,3), {Boundary => Reflect}

=for options

 Boundary - controls what values are assumed for the image when kernel
            crosses its edge:
 	    => Default   - periodic boundary conditions
                           (i.e. wrap around axis)
 	    => Reflect   - reflect at boundary
 	    => Truncate  - truncate at boundary
 	    => Replicate - repeat boundary pixel values

=cut
EOD
       BadDoc =>
'Unlike the FFT routines, conv2d is able to process bad values.',
       HandleBad => 1,
        Pars => 'a(m,n); kern(p,q); [o]b(m,n); indx [t]mapi(isize=CALC($SIZE(p) + $SIZE(m))); indx [t]mapj(jsize=CALC($SIZE(q) + $SIZE(n)))',
        OtherPars => 'int opt;',
        PMCode => '
sub PDL::conv2d {
   my $opt; $opt = pop @_ if ref($_[$#_]) eq \'HASH\';
   die \'Usage: conv2d( a(m,n), kern(p,q), [o]b(m,n), {Options} )\'
      if $#_<1 || $#_>2;
   my($x,$kern) = @_;
   my $c = $#_ == 2 ? $_[2] : $x->nullcreate;
   PDL::_conv2d_int($x,$kern,$c,
	(!($opt && exists $$opt{Boundary}))?0:$boundary2value{$$opt{Boundary}}
	);
   return $c;
}
',
        GenericTypes => $A,
        Code =>
           init_vars() .
           init_map("i") .
           init_map("j") .
pp_line_numbers(__LINE__, <<'EOF'),
broadcastloop %{
loop(n,m) %{
  PDL_IF_BAD(int flag = 0;,)
  PDL_CLDouble tmp = 0;
  loop (q) %{
     j2 = mapj[n-q];
     if (j2 >= 0) {
        loop (p) %{
           i2 = mapi[m-p];
           if (i2 >= 0) {
              PDL_IF_BAD(if ( $ISGOOD(a(m=>i2,n=>j2)) && $ISGOOD(kern()) ),) {
                 tmp += $a(m=>i2,n=>j2) * $kern();
                 PDL_IF_BAD(flag = 1;,)
              } /* if: good */
           } /* if: i2 >= 0 */
        %}
     } /* if: j2 >= 0 */
  %}
  PDL_IF_BAD(if ( !flag ) { $SETBAD(b()); }
  else,)                  { $b() = tmp; }
%}
%}
EOF
); # pp_def: conv2d

pp_def('med2d', Doc=> <<'EOD',
=for ref

2D median-convolution of an array with a kernel (smoothing)

Note: only points in the kernel E0 are included in the median, other
points are weighted by the kernel value (medianing lots of zeroes
is rather pointless)

=for usage

 $new = med2d $old, $kernel, {OPTIONS}

=for example

 $smoothed = med2d $image, ones(3,3), {Boundary => Reflect}

=for options

 Boundary - controls what values are assumed for the image when kernel
            crosses its edge:
 	    => Default   - periodic boundary conditions (i.e. wrap around axis)
 	    => Reflect   - reflect at boundary
 	    => Truncate  - truncate at boundary
 	    => Replicate - repeat boundary pixel values

=cut
EOD
       BadDoc =>
'Bad values are ignored in the calculation. If all elements within the
kernel are bad, the output is set bad.',
       HandleBad => 1,
        Pars => 'a(m,n); kern(p,q); [o]b(m,n); double+ [t]tmp(pq=CALC($SIZE(p)*$SIZE(q))); indx [t]mapi(isize=CALC($SIZE(p) + $SIZE(m))); indx [t]mapj(jsize=CALC($SIZE(q) + $SIZE(n)))',
        OtherPars => 'int opt;',
        PMCode => pp_line_numbers(__LINE__, <<'EOF'),
sub PDL::med2d {
   my $opt; $opt = pop @_ if ref($_[$#_]) eq 'HASH';
   die 'Usage: med2d( a(m,n), kern(p,q), [o]b(m,n), {Options} )'
      if $#_<1 || $#_>2;
   my($x,$kern) = @_;
   croak "med2d: kernel must contain some positive elements.\n"
       if all( $kern <= 0 );
   my $c = $#_ == 2 ? $_[2] : $x->nullcreate;
   PDL::_med2d_int($x,$kern,$c,
	(!($opt && exists $$opt{Boundary}))?0:$boundary2value{$$opt{Boundary}}
   );
   return $c;
}
EOF
        CHeader => 'PDL_TYPELIST_FLOATREAL(PDL_QSORT)',
        Code =>
           init_vars() .
           init_map("i") .
           init_map("j") .
pp_line_numbers(__LINE__, <<'EOF'),
broadcastloop %{
loop (n,m) %{
  PDL_Indx count = 0;
  PDL_IF_BAD(int flag = 0;,)
  loop (q) %{
    j2 = mapj[n-q];
    if (j2 >= 0)
      loop (p) %{
        i2 = mapi[m-p];
        if (i2 >= 0) {
          PDL_LDouble kk = $kern(), aa = $a(m=>i2,n=>j2);
          PDL_IF_BAD(if ( $ISGOODVAR(kk,kern) && $ISGOODVAR(aa,a) ),) {
            PDL_IF_BAD(flag = 1;,)
            if ( kk > 0 )
              $tmp(pq=>count++) = aa * kk;
          }
        } /* if: i2 >= 0 */
      %}
  %}
  PDL_IF_BAD(if ( flag == 0 ) {
    $SETBAD(b());
  } else,) {
    qsort_$PPSYM(tmp)( $P(tmp), 0, count-1 );
    $b() = $tmp(pq=>(count-1)/2);
  }
%}
%}
EOF
); # pp_def: med2d

pp_def('med2df', Doc=> <<'EOD',

=for ref

2D median-convolution of an array in a pxq window (smoothing)

Note: this routine does the median over all points in a rectangular
      window and is not quite as flexible as C in this regard
      but slightly faster instead

=for usage

 $new = med2df $old, $xwidth, $ywidth, {OPTIONS}

=for example

 $smoothed = med2df $image, 3, 3, {Boundary => Reflect}

=for options

 Boundary - controls what values are assumed for the image when kernel
            crosses its edge:
 	    => Default   - periodic boundary conditions (i.e. wrap around axis)
 	    => Reflect   - reflect at boundary
 	    => Truncate  - truncate at boundary
 	    => Replicate - repeat boundary pixel values

=cut

EOD
        Pars => 'a(m,n); [o]b(m,n); indx [t]mapi(isize=CALC($SIZE(p) + $SIZE(m))); indx [t]mapj(jsize=CALC($SIZE(q) + $SIZE(n)))',
        OtherPars => 'IV p_size=>p; IV q_size=>q; int opt;',
        PMCode => <<'EOF',
sub PDL::med2df {
   my $opt; $opt = pop @_ if ref($_[$#_]) eq 'HASH';
   die 'Usage: med2df( a(m,n), [o]b(m,n), p, q, {Options} )'
      if $#_<2 || $#_>3;
   my($x,$p,$q) = @_;
   croak "med2df: kernel must contain some positive elements.\n"
       if $p == 0 && $q == 0;
   my $c = $#_ == 3 ? $_[3] : $x->nullcreate;
   &PDL::_med2df_int($x,$c,$p,$q,
	(!($opt && exists $$opt{Boundary}))?0:$boundary2value{$$opt{Boundary}}
   );
   return $c;
}
EOF
        Code =>
           init_vars() .
           init_map("i") .
           init_map("j") .
'
  $GENERIC() tmp[p_size*q_size];
  broadcastloop %{
  loop (n,m) %{
     int count = 0;
     loop (q) %{
        if ((j2 = mapj[n-q]) < 0) continue;
        loop (p) %{
           if ((i2 = mapi[m-p]) < 0) continue;
           tmp[count++] = $a(m=>i2,n=>j2);
        %}
     %}
     $b() = quick_select_$PPSYM() (tmp, count );
  %}
  %}
',

); # pp_def: med2df

pp_addhdr(<<'EOH');
#define EZ(x) ez ? 0 : (x)
EOH
pp_def('box2d',
       Pars => 'a(n,m); [o] b(n,m)',
       OtherPars => 'int wx; int wy; int edgezero',
       GenericTypes => $A,
       Code => '
                register int nx = 0.5*$COMP(wx);
                register int ny = 0.5*$COMP(wy);
		register int ez = $COMP(edgezero);
                PDL_CLDouble div, sum, lsum;
                int ind1,ind2,first;

		div = 1/((2.0*nx+1)*(2.0*ny+1));

		broadcastloop %{
		  first = 1;
		  loop (m,n=:nx) %{
		    ind1 = $SIZE(n)-1-n;  /* rightmost strip */
		    $b() = EZ($a());
		    $b(n=>ind1) = EZ($a(n=>ind1));
                  %}
		  loop (n,m=:ny) %{
		    ind1 = $SIZE(m)-1-m;  /* topmost strip */
		    $b() = EZ($a());
		    $b(m=>ind1) = EZ($a(m=>ind1));
		  %}
		  loop (m=ny:-ny) %{
		    if (first) {
		       first = 0;
		       lsum = 0;
		       PDL_Indx m_outer = m;
		       loop (m=m_outer-ny:m_outer+ny+1,n=:2*nx+1) %{
			 lsum += $a(); /* initialize sum */
		       %}
		    } else {
			ind1 = m-ny-1;
			ind2 = m+ny;
			loop (n=:2*nx+1) %{
			  lsum -= $a(m=>ind1); /* remove top pixels */
			  lsum += $a(m=>ind2); /* add bottom pixels */
			%}
		    }
		    sum = lsum;
		    $b(n=>nx) = div*sum;     /* and assign */
		    loop (n=nx+1:-nx) %{ /* loop along line */
		      ind1 = n-nx-1;
		      ind2 = n+nx;
		      PDL_Indx m_outer = m;
		      loop (m=m_outer-ny:m_outer+ny+1) %{
			sum -= $a(n=>ind1); /* remove leftmost data */
			sum += $a(n=>ind2); /* and add rightmost */
		      %}
		      $b() = div*sum;         /* and assign */
		    %}
		  %}
		%}',
       Doc => << 'EOD',
=for ref

fast 2D boxcar average

=for example

  $smoothim = $im->box2d($wx,$wy,$edgezero=1);

The edgezero argument controls if edge is set to zero (edgezero=1)
or just keeps the original (unfiltered) values.

C should be updated to support similar edge options
as C and C etc.

Boxcar averaging is a pretty crude way of filtering. For serious stuff
better filters are around (e.g., use L with the appropriate
kernel). On the other hand it is fast and computational cost grows only
approximately linearly with window size.

=cut

EOD
); # pp_def box2d

=head2 patch2d

=cut

pp_def('patch2d',
       Doc=><<'EOD',

=for ref

patch bad pixels out of 2D images using a mask

C<$bad> is a 2D mask array where 1=bad pixel 0=good pixel.
Pixels are replaced by the average of their non-bad neighbours;
if all neighbours are bad, the original data value is
copied across.

=cut

EOD
       BadDoc =>
'This routine does not handle bad values - use L instead',
       HandleBad => 0,
        Pars => 'a(m,n); int bad(m,n); [o]b(m,n);',
        GenericTypes => $A,
        Code =>
        'PDL_Indx i1,j1, i2,j2, norm;
         PDL_CLDouble tmp;
      broadcastloop %{
         loop (n,m) %{
           $b() = $a();
           if (!$bad()) continue;
           tmp = 0; norm=0;
           for(j1=-1; j1<=1; j1++) {
              j2 = n+j1;
              if ( j2<0 || j2>=$SIZE(n) ) continue;
              for(i1=-1; i1<=1; i1++) {
                 /* ignore central pixel, which we know is bad */
                 if ( i1==0 && j1==0 ) continue;
                 i2 = m+i1;
                 if ( i2<0 || i2>=$SIZE(m) || $bad(m=>i2,n=>j2) ) continue;
                 tmp += $a(m=>i2,n=>j2);
                 norm++;
              } /* for: i1 */
           } /* for: j1 */
           if (norm>0) {  /* Patch */
              $b() = tmp/norm;
           }
         %}
      %} /* broadcastloop */
         ', # Code
);

pp_def('patchbad2d',
       Doc=><<'EOD',
=for ref

patch bad pixels out of 2D images containing bad values

Pixels are replaced by the average of their non-bad neighbours;
if all neighbours are bad, the output is set bad.
If the input ndarray contains I bad values, then a straight copy
is performed (see L).
EOD
       BadDoc =>
'patchbad2d handles bad values. The output ndarray I contain
bad values, depending on the pattern of bad values in the input ndarray.',
       HandleBad => 1,
        Pars => 'a(m,n); [o]b(m,n);',
        GenericTypes => $A,
        Code => pp_line_numbers(__LINE__, <<'EOF'),
PDL_IF_BAD(int flag = 0;,)
broadcastloop %{
   loop(n,m) %{
      $GENERIC(a) a_val = $a();
      PDL_IF_BAD(if ( !$ISGOODVAR(a_val,a) ) {
         PDL_CLDouble tmp = 0; PDL_Indx norm=0; int i1; int j1;
         for(j1=-1; j1<=1; j1++) {
            PDL_Indx j2 = n+j1;
            if ( j2<0 || j2>=$SIZE(n) ) continue;
            for(i1=-1; i1<=1; i1++) {
               /* ignore central pixel, which we know is bad */
               if ( i1==0 && j1==0 ) continue;
               PDL_Indx i2 = m+i1;
               if ( i2<0 || i2>=$SIZE(m) ) continue;
               a_val = $a(m=>i2,n=>j2);
               if ( $ISBADVAR(a_val,a) ) continue;
               tmp += a_val;
               norm++;
            } /* for: i1 */
         } /* for: j1 */
         /* Patch */
         if (norm>0) {
            $b() = tmp/norm;
         } else {
            $SETBAD(b());
            flag = 1;
         }
      } else,) {
         $b() = a_val;
      } /* if: ISGOODVAR() */
   %}
%} /* broadcastloop */
/* handle bad flag */
PDL_IF_BAD(if ( flag ) $PDLSTATESETBAD(b);,)
EOF
);

pp_def('max2d_ind',
      Doc=><<'EOD',

=for ref

Return value/position of maximum value in 2D image

Contributed by Tim Jenness

=cut

EOD

      BadDoc=><<'EOD',

Bad values are excluded from the search. If all pixels
are bad then the output is set bad.

EOD

       HandleBad => 1,
        Pars => 'a(m,n); [o]val(); int [o]x(); int[o]y();',
        Code => '
        PDL_LDouble cur = 0;
        PDL_Indx curind1 = PDL_IF_BAD(-1,0), curind2 = PDL_IF_BAD(-1,0);
        loop(m,n) %{
          if (!(PDL_IF_BAD($ISGOOD(a()) &&,) ( (!n && !m) || $a() > cur || isnan(cur)))) continue;
          cur = $a(); curind1 = m; curind2 = n;
        %}
        PDL_IF_BAD(if ( curind1 < 0 ) {
          $SETBAD(val());
          $SETBAD(x());
          $SETBAD(y());
        } else,) {
          $val() = cur;
          $x()   = curind1;
          $y()   = curind2;
	}
        ');

pp_def('centroid2d',
        Doc=><<'EOD',
=for ref

Refine a list of object positions in 2D image by centroiding in a box

C<$box> is the full-width of the box, i.e. the window
is C<+/- $box/2>.

=cut
EOD
        BadDoc=><<'EOD',
Bad pixels are excluded from the centroid calculation. If all elements are
bad (or the pixel sum is 0 - but why would you be centroiding
something with negatives in...) then the output values are set bad.
EOD
       HandleBad => 1,
	Pars => 'im(m,n); x(); y(); box(); [o]xcen(); [o]ycen();',
	Code => '
PDL_LDouble sum = 0, sumx = 0, sumy = 0;
PDL_Indx i1 = $x() - $box()/2;
PDL_Indx i2 = $x() + $box()/2 + 1;
PDL_Indx j1 = $y() - $box()/2;
PDL_Indx j2 = $y() + $box()/2 + 1;
loop (n=j1:j2,m=i1:i2) %{
   PDL_LDouble data = $im();
   PDL_IF_BAD(if ( $ISBADVAR(data,im) ) continue;,)
   sum += data;
   sumx += data*m;
   sumy += data*n;
%}
/*
 * if sum == 0 then we will flag as bad -- although it could just mean that
 * there is negative values in the dataset.
 * - should use a better check than != 0.0  ...
 */
PDL_IF_BAD(if ( sum == 0.0 ) {
   $SETBAD(xcen());
   $SETBAD(ycen());
} else,) {
   $xcen() = sumx/sum;
   $ycen() = sumy/sum;
}
'
);

pp_add_exported('', 'crop');
pp_addpm(<<'EOPM');
=head2 crop

=for ref

Return bounding box of given mask in an C ndarray, so it can broadcast.
Use other operations (such as L, or
L with a colour vector) to create a mask suitable
for your application.

=for example

  $x1x2y1y2 = crop($image);

=cut

*crop = \&PDL::crop;
sub PDL::crop {
  my ($mask) = @_;
  $mask->xchg(0,1)->orover->_which_int(my $out = null, null);
  $out->badflag(1); $out->badvalue(-1);
  my ($x1, $x2) = $out->minmaximum;
  $mask->orover->_which_int($out = null, null);
  $out->badflag(1); $out->badvalue(-1);
  my ($y1, $y2) = $out->minmaximum;
  $x1->cat($x2, $y1, $y2)->mv(-1,0);
}
EOPM

pp_add_exported('', 'cc8compt','cc4compt');
pp_addpm(<<'EOPM');

=head2 cc8compt

=for ref

Connected 8-component labeling of a binary image.

Connected 8-component labeling of 0,1 image - i.e. find separate
segmented objects and fill object pixels with object number.
8-component labeling includes all neighboring pixels.
This is just a front-end to ccNcompt.  See also L.

=for example

 $segmented = cc8compt( $image > $threshold );

=head2 cc4compt

=for ref

Connected 4-component labeling of a binary image.

Connected 4-component labeling of 0,1 image - i.e. find separate
segmented objects and fill object pixels with object number.
4-component labling does not include the diagonal neighbors.
This is just a front-end to ccNcompt.  See also L.

=for example

 $segmented = cc4compt( $image > $threshold );

=cut

sub PDL::cc8compt{
return ccNcompt(shift,8);
}
*cc8compt = \&PDL::cc8compt;

sub PDL::cc4compt{
return ccNcompt(shift,4);
}
*cc4compt = \&PDL::cc4compt;

EOPM

pp_def('ccNcompt',Doc=>'

=for ref

Connected component labeling of a binary image.

Connected component labeling of 0,1 image - i.e. find separate
segmented objects and fill object pixels with object number.
See also L and L.

The connectivity parameter must be 4 or 8.

=for example

 $segmented = ccNcompt( $image > $threshold, 4);

 $segmented2 = ccNcompt( $image > $threshold, 8);

where the second parameter specifies the connectivity (4 or 8) of the labeling.

=cut

',
       HandleBad => 0, # a marker
        Pars => 'a(m,n); int+ [o]b(m,n);',
        OtherPars => 'int con',
        CHeader => qq{void AddEquiv(PDL_Long* equiv, PDL_Long i, PDL_Long j);\n},
        Code => '

      PDL_Long i,k;
      PDL_Long newlabel;
      PDL_Long neighbour[4];
      PDL_Long nfound;
      PDL_Long pass,count,next,this;
      PDL_Long *equiv = NULL;
      PDL_Long i1,j1,i2;

      if ($COMP(con)!=4 && $COMP(con)!=8)
         $CROAK("In ccNcompt, connectivity must be 4 or 8, you gave %d",$COMP(con));
      loop(n,m) %{ /* Copy */
         $b() = $a();
      %}

      /* 1st pass counts max possible compts, 2nd records equivalences */

      for (pass = 0; pass<2; pass++) {

      if (pass==1) {
         equiv = (PDL_Long*) malloc((newlabel+1)*sizeof(PDL_Long));
         if (equiv==(PDL_Long*)0)
            $CROAK("Out of memory");
         for(i=0;i<=newlabel;i++)
             equiv[i]=i;
      }

      newlabel = 1; /* Running label */

      loop (n,m) %{ /* Loop over image pixels */
            nfound = 0; /* Number of neighbour >0 */

            i1 = m-1; j1 = n-1; i2 = m+1; /*West x, North y, East x*/

            if ($b() > 0) { /* Check 4 neighbour already seen */

               if (m>0 && $b(m=>i1)>0) /*West*/
                   neighbour[nfound++] = $b(m=>i1); /* Store label of it */
               if (n>0 && $b(n=>j1)>0) /*North*/
                   neighbour[nfound++] = $b(n=>j1);
               if (n>0 && m>0  && $b(m=>i1, n=>j1)>0 && $COMP(con) == 8) /*North-West*/
                   neighbour[nfound++] = $b(m=>i1, n=>j1);
               if (n>0 && m<($SIZE(m)-1) && $b(m=>i2, n=>j1)>0 && $COMP(con) == 8) /*North-East*/
                   neighbour[nfound++] = $b(m=>i2, n=>j1);

               if (nfound==0)  { /* Assign new label */
                  $b() = newlabel++;
               }
               else {
                  $b() =  neighbour[0];
                  if (nfound>1 && pass == 1) {  /* Assign equivalents */
                      for(k=1; k
                         AddEquiv( equiv, (PDL_Long)$b(),
                            neighbour[k] );
                  }
               }
            }

            else {  /* No label */

               $b() = 0;
            }

      %} /* End of image loop */

      } /* Passes */

      /* Replace each cycle by single label */

       count = 0;
       for (i = 1; i <= newlabel; i++)
         if ( i <= equiv[i] ) {
             count++;
             this = i;
             while ( equiv[this] != i ) {
               next = equiv[this];
               equiv[this] = count;
               this = next;
             }
          equiv[this] = count;
         }


      /* Now remove equivalences */

      loop (n,m) %{ /* Loop over image pixels */
           $b()   = equiv[ (PDL_Long) $b()  ] ;
      %}

      free(equiv); /* Tidy */
');

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

=for ref

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

This function works inplace, i.e. modifies C.

=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. Often, it includes vertices
and edge points. Set the C 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 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');

=head2 pnpoly

=for ref

'points in a polygon' selection from a 2-D ndarray

=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 {$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
L
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 >.

  $px = pdl( 3, 20, 34 );
  $py = pdl( 3, 20,  3 );
  $x = $img->xvals;      # get x pixel coords
  $y = $img->yvals;      # get y pixel coords

  # $tri is 0 everywhere except for points in polygon interior
  $tri = pnpoly($x,$y,$px,$py);

=cut

# From: http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
#
# Fixes needed to pnpoly code:
#
# Use topdl() to ensure ndarray args
#
# Add POD docs for usage
#
# Calculate first term in & expression to use to fix divide-by-zero when
#   the test point is in line with a vertical edge of the polygon.
#   By adding the value of $mask we prevent a divide-by-zero since the &
#   operator does not "short circuit".

sub PDL::pnpoly {
	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_add_exported('', 'polyfillv');
pp_addpm(<<'EOPM');

=head2 polyfillv

=for ref

return the (dataflowed) 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. Often, it includes vertices
and edge points. Set the C 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 algorithm

=for example

  # 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 $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';
	}

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

EOPM

pp_addhdr(<<'EOF');
int getnewsize(PDL_Indx cols, PDL_Indx rows, float fangle, PDL_Indx *newcols, PDL_Indx *newrows);
typedef unsigned char imT;    /* image type */
int rotate(imT *im, imT *out, int cols, int rows, int nc, int nr,
           float fangle, imT bgval, int antialias);
EOF
pp_add_exported('','rotnewsz');
pp_addxs('

void
rotnewsz(m,n,angle)
	int m
	int n
	float angle
	PPCODE:
	PDL_Indx newcols, newrows;

	if (getnewsize(m,n,angle,&newcols,&newrows) != 0)
		croak("wrong angle (should be between -90 and +90)");
	EXTEND(sp,2);
	PUSHs(sv_2mortal(newSVnv(newcols)));
	PUSHs(sv_2mortal(newSVnv(newrows)));
');

pp_def('rot2d',
       HandleBad => 0,
	Pars => 'im(m,n); float angle(); bg(); int aa(); [o] om(p,q)',
	GenericTypes => $A,
	Code => 'int ierr;
		 if ((ierr = rotate($P(im),$P(om),$SIZE(m),$SIZE(n),$SIZE(p),
			$SIZE(q),$angle(),$bg(),$aa())) != 0) {
			if (ierr == -1)
				$CROAK("error during rotate, wrong angle");
			else
				$CROAK("wrong output dims, did you set them?");
		}',
	RedoDimsCode => '
			if (getnewsize($SIZE(m),$SIZE(n),
				$angle(), &$SIZE(p),
				&$SIZE(q)) != 0)
			   $CROAK("error during rotate, wrong angle");
			/* printf("o: %d, p: %d\n",ncols,nrows); */
			',
	GenericTypes => ['B'],
	Doc => << 'EOD',

=for ref

rotate an image by given C

=for example

  # rotate by 10.5 degrees with antialiasing, set missing values to 7
  $rot = $im->rot2d(10.5,7,1);

This function rotates an image through an C between -90 and + 90
degrees. Uses/doesn't use antialiasing depending on the C flag.
Pixels outside the rotated image are set to C.

Code modified from pnmrotate (Copyright Jef Poskanzer) with an algorithm based
on "A Fast Algorithm for General  Raster  Rotation"  by  Alan Paeth,
Graphics Interface '86, pp. 77-81.

Use the C function to find out about the dimension of the
newly created image

  ($newcols,$newrows) = rotnewsz $oldn, $oldm, $angle;

L offers a more general interface to
distortions, including rotation, with various types of sampling; but
rot2d is faster.

=cut

EOD
);

pp_def('bilin2d',
       HandleBad => 0,
    Pars => 'Int(n,m); [io] O(q,p)',
    Doc=><<'EOD',
=for ref

Bilinearly maps the first ndarray in the second. The
interpolated values are actually added to the second
ndarray which is supposed to be larger than the first one.
EOD
    GenericTypes => $A,
    RedoDimsCode =>'
     if ($SIZE(q)<$SIZE(n) || $SIZE(p)<$SIZE(m))
       $CROAK("the second matrix must be greater than first!");
    ',
    Code =>'
  PDL_CLDouble dx = ((PDL_CLDouble) ($SIZE(n)-1)) / ($SIZE(q)-1);
  PDL_CLDouble dy = ((PDL_CLDouble) ($SIZE(m)-1)) / ($SIZE(p)-1);
  broadcastloop %{
    PDL_CLDouble x=0;
    loop (q) %{
      PDL_CLDouble y = 0;
      loop (p) %{
        PDL_Indx ii = floor(x);
        if (ii>=($SIZE(n)-1)) ii = $SIZE(n)-2;
        PDL_Indx jj = floor(y);
        if (jj>=($SIZE(m)-1)) jj = $SIZE(m)-2;
        PDL_Indx ii1 = ii+1;
        PDL_Indx jj1 = jj+1;
        PDL_CLDouble y1 = $Int(n=>ii,m=>jj);
        PDL_CLDouble y2 = $Int(n=>ii1,m=>jj);
        PDL_CLDouble y3 = $Int(n=>ii1,m=>jj1);
        PDL_CLDouble y4 = $Int(n=>ii,m=>jj1);
        PDL_CLDouble t = x-ii;
        PDL_CLDouble u = y-jj;
        $O() += (1-t)*(1-u)*y1 + t*(1-u)*y2 + t*u*y3 + (1-t)*u*y4;
        y+=dy;
      %}
      x+=dx;
    %}
  %}
');

pp_def('rescale2d',
       HandleBad => 0,
    Pars => 'Int(m,n); [io] O(p,q)',
    Doc=><<'EOD',
=for ref

The first ndarray is rescaled to the dimensions of the second
(expanding or meaning values as needed) and then added to it in place.
Nothing useful is returned.

If you want photometric accuracy or automatic FITS header metadata
tracking, consider using L
instead: it does these things, at some speed penalty compared to
rescale2d.
EOD
    GenericTypes => $A,
    Code =>'
if($SIZE(p) >= $SIZE(m) && $SIZE(q) >= $SIZE(n)) {
  PDL_LDouble kx = ((PDL_LDouble) $SIZE(p)) / $SIZE(m);
  PDL_LDouble ky = ((PDL_LDouble) $SIZE(q)) / $SIZE(n);
  broadcastloop %{
    PDL_Indx lx = 0;
    loop (m) %{
      PDL_Indx ly = 0, cx = rint((m+1)*kx);
      loop (n) %{
	PDL_Indx cy = rint((n+1)*ky);
	loop (p=lx:cx,q=ly:cy) %{
	  $O() += $Int();
	%}
	ly = cy + 1;
      %}
      lx = cx + 1;
    %}
  %}
}
else if($SIZE(p) < $SIZE(m) && $SIZE(q) < $SIZE(n)) {
  PDL_LDouble kx = ((PDL_LDouble) $SIZE(m)) / $SIZE(p);
  PDL_LDouble ky = ((PDL_LDouble) $SIZE(n)) / $SIZE(q);
  broadcastloop %{
    PDL_Indx lx = 0;
    loop (p) %{
      PDL_Indx ly = 0, cx = rint((p+1)*kx);
      loop (q) %{
	PDL_Indx cy = rint((q+1)*ky);
	PDL_LDouble temp = 0.0;
	PDL_Indx num = 0;
	loop (m=lx:cx,n=ly:cy) %{
	  temp += $Int();
	  num++;
	%}
	$O() += temp/num;
	ly = cy + 1;
      %}
      lx = cx + 1;
    %}
  %}
}
else if($SIZE(p) >= $SIZE(m) && $SIZE(q) < $SIZE(n)) {
  PDL_LDouble kx = ((PDL_LDouble) $SIZE(p)) / $SIZE(m);
  PDL_LDouble ky = ((PDL_LDouble) $SIZE(n)) / $SIZE(q);
  broadcastloop %{
    PDL_Indx lx = 0;
    loop (m) %{
      PDL_Indx ly = 0, cx = rint((m+1)*kx);
      loop (q) %{
	PDL_Indx cy = rint((q+1)*ky);
	PDL_LDouble temp = 0.0;
	PDL_Indx num = 0;
	loop (n=ly:cy) %{
	  temp += $Int();
	  num++;
	%}
	loop (p=lx:cx) %{
	  $O() += temp/num;
	%}
	ly = cy + 1;
      %}
      lx = cx + 1;
    %}
  %}
}
else if($SIZE(p) < $SIZE(m) && $SIZE(q) >= $SIZE(n)) {
  PDL_LDouble kx = ((PDL_LDouble) $SIZE(m)) / $SIZE(p);
  PDL_LDouble ky = ((PDL_LDouble) $SIZE(q)) / $SIZE(n);
  broadcastloop %{
    PDL_Indx lx = 0;
    loop (p) %{
      PDL_Indx ly = 0, cx = rint((p+1)*kx);
      loop (n) %{
	PDL_Indx cy = rint((n+1)*ky);
	PDL_CLDouble temp = 0.0;
	PDL_Indx num = 0;
	loop (m=lx:cx) %{
	  temp += $Int();
	  num++;
	%}
	loop (q=ly:cy) %{
	  $O() += temp/num;
	%}
	ly = cy + 1;
      %}
      lx = cx + 1;
    %}
  %}
}
else $CROAK("I am not supposed to be here, please report the bug to ");
  ');

# functions to make handling 2D polynomial mappings a bit easier
#
pp_add_exported('', 'fitwarp2d applywarp2d');

pp_addpm( '

=head2 fitwarp2d

=for ref

Find the best-fit 2D polynomial to describe
a coordinate transformation.

=for usage

  ( $px, $py ) = fitwarp2d( $x, $y, $u, $v, $nf, { options } )

Given a set of points in the output plane (C<$u,$v>), find
the best-fit (using singular-value decomposition) 2D polynomial
to describe the mapping back to the image plane (C<$x,$y>).
The order of the fit is controlled by the C<$nf> parameter
(the maximum power of the polynomial is C<$nf - 1>), and you
can restrict the terms to fit using the C option.

C<$px> and C<$py> are C by C element ndarrays which describe
a polynomial mapping (of order C)
from the I C<(u,v)> image to the I C<(x,y)> image:

  x = sum(j=0,np-1) sum(i=0,np-1) px(i,j) * u^i * v^j
  y = sum(j=0,np-1) sum(i=0,np-1) py(i,j) * u^i * v^j

The transformation is returned for the reverse direction (ie
output to input image) since that is what is required by the
L routine.  The L
routine can be used to convert a set of C<$u,$v> points given
C<$px> and C<$py>.

Options:

=for options

  FIT     - which terms to fit? default ones(byte,$nf,$nf)

=begin comment

old option, caused trouble
  THRESH  - in svd, remove terms smaller than THRESH * max value
            default is 1.0e-5

=end comment

=over 4

=item FIT

C allows you to restrict which terms of the polynomial to fit:
only those terms for which the FIT ndarray evaluates to true will be
evaluated.  If a 2D ndarray is sent in, then it is
used for the x and y polynomials; otherwise
C<< $fit->slice(":,:,(0)") >> will be used for C<$px> and
C<< $fit->slice(":,:,(1)") >> will be used for C<$py>.

=begin comment

=item THRESH

Remove all singular values whose value is less than C
times the largest singular value.

=end comment

=back

The number of points must be at least equal to the number of
terms to fit (C<$nf*$nf> points for the default value of C).

=for example

  # points in original image
  $x = pdl( 0,   0, 100, 100 );
  $y = pdl( 0, 100, 100,   0 );
  # get warped to these positions
  $u = pdl( 10, 10, 90, 90 );
  $v = pdl( 10, 90, 90, 10 );
  #
  # shift of origin + scale x/y axis only
  $fit = byte( [ [1,1], [0,0] ], [ [1,0], [1,0] ] );
  ( $px, $py ) = fitwarp2d( $x, $y, $u, $v, 2, { FIT => $fit } );
  print "px = ${px}py = $py";
  px =
  [
   [-12.5  1.25]
   [    0     0]
  ]
  py =
  [
   [-12.5     0]
   [ 1.25     0]
  ]
  #
  # Compared to allowing all 4 terms
  ( $px, $py ) = fitwarp2d( $x, $y, $u, $v, 2 );
  print "px = ${px}py = $py";
  px =
  [
   [         -12.5           1.25]
   [  1.110223e-16 -1.1275703e-17]
  ]
  py =
  [
   [         -12.5  1.6653345e-16]
   [          1.25 -5.8546917e-18]
  ]

  # A higher-degree polynomial should not affect the answer much, but
  # will require more control points

  $x = $x->glue(0,pdl(50,12.5, 37.5, 12.5, 37.5));
  $y = $y->glue(0,pdl(50,12.5, 37.5, 37.5, 12.5));
  $u = $u->glue(0,pdl(73,20,40,20,40));
  $v = $v->glue(0,pdl(29,20,40,40,20));
  ( $px3, $py3 ) = fitwarp2d( $x, $y, $u, $v, 3 );
  print "px3 =${px3}py3 =$py3";
  px3 =
  [
   [-6.4981162e+08       71034917     -726498.95]
   [      49902244     -5415096.7      55945.388]
   [    -807778.46      88457.191     -902.51612]
  ]
  py3 =
  [
   [-6.2732159e+08       68576392     -701354.77]
   [      48175125     -5227679.8      54009.114]
   [    -779821.18      85395.681     -871.27997]
  ]

  #This illustrates an important point about singular value
  #decompositions that are used in fitwarp2d: like all SVDs, the
  #rotation matrices are not unique, and so the $px and $py returned
  #by fitwarp2d are not guaranteed to be the "simplest" solution.
  #They do still work, though:

  ($x3,$y3) = applywarp2d($px3,$py3,$u,$v);
  print approx $x3,$x,1e-4;
  [1 1 1 1 1 1 1 1 1]
  print approx $y3,$y;
  [1 1 1 1 1 1 1 1 1]

=head2 applywarp2d

=for ref

Transform a set of points using a 2-D polynomial mapping

=for usage

  ( $x, $y ) = applywarp2d( $px, $py, $u, $v )

Convert a set of points (stored in 1D ndarrays C<$u,$v>)
to C<$x,$y> using the 2-D polynomial with coefficients stored in C<$px>
and C<$py>.  See L
for more information on the format of C<$px> and C<$py>.

=cut

# use SVD to fit data. Assuming no errors.

=pod

=begin comment

Some explanation of the following three subroutines (_svd, _mkbasis,
and fitwarp2d): See Wolberg 1990 (full ref elsewhere in this
documentation), Chapter 3.6 "Polynomial Transformations".  The idea is
that, given a set of control points in the input and output images
denoted by coordinates (x,y) and (u,v), we want to create a polynomial
transformation that relates u to linear combinations of powers of x
and y, and another that relates v to powers of x and y.

The transformations used here and by Wolberg differ slightly, but the
basic idea is something like this.  For each u and each v, define a
transform:

u = (sum over j) (sum over i) a_{ij} x**i * y**j , (eqn 1)
v = (sum over j) (sum over i) b_{ij} x**i * y**j . (eqn 2)

We can write this in matrix notation.  Given that there are M control
points, U is a column vector with M rows.  A and B are vectors containing
the N coefficients (related to the degree of the polynomial fit).  W
is an MxN matrix of the basis row-vectors (the x**i y**j).

The matrix equations we are trying to solve are
U = W A , (eqn 3)
V = W B . (eqn 4)

We need to find the A and B column vectors, those are the coefficients
of the polynomial terms in W.  W is not square, so it has no inverse.
But is has a pseudo-inverse W^+ that is NxM.  We are going to use that
pseudo-inverse to isolate A and B, like so:

W^+ U = W^+ W A = A (eqn 5)
W^+ V = W^+ W B = B (eqn 6)

We are going to get W^+ by performing a singular value decomposition
of W (to use some of the variable names below):

W = $svd_u x SIGMA x $svd_v->transpose (eqn 7)
W^+ = $svd_v x SIGMA^+ x $svd_u->transpose . (eqn 8)

Here SIGMA is a square diagonal matrix that contains the singular
values of W that are in the variable $svd_w.  SIGMA^+ is the
pseudo-inverse of SIGMA, which is calculated by replacing the non-zero
singular values with their reciprocals, and then taking the transpose
of the matrix (which is a no-op since the matrix is square and
diagonal).

So the code below does this:

_mkbasis computes the matrix W, given control coordinates u and v and
the maximum degree of the polynomial (and the terms to use).

_svd takes the svd of W, computes the pseudo-inverse of W, and then
multiplies that with the U vector (there called $y). The output of
_svd is the A or B vector in eqns 5 & 6 above. Rarely is the matrix
multiplication explicit, unfortunately, so I have added EXPLANATIONs
below.

=end comment

=cut

sub _svd ($$) {
    my $basis  = shift;
    my $y      = shift;
#    my $thresh = shift;

    # if we had errors for these points, would normalise the
    # basis functions, and the output array, by these errors here

    # perform the SVD
    my ( $svd_u, $svd_w, $svd_v ) = svd( $basis );

    # DAL, 09/2017: the reason for removing ANY singular values, much less
    #the smallest ones, is not clear. I commented the line below since this
    #actually removes the LARGEST values in SIGMA^+.
    # $svd_w *= ( $svd_w >= ($svd_w->max * $thresh ) );
    # The line below would instead remove the SMALLEST values in SIGMA^+, but I can see no reason to include it either.
    # $svd_w *= ( $svd_w <= ($svd_w->min / $thresh ) );

    # perform the back substitution
    # EXPLANATION: same thing as $svd_u->transpose x $y->transpose.
    my $tmp = $y x $svd_u;

    #EXPLANATION: the division by (the non-zero elements of) $svd_w (the singular values) is
    #equivalent to $sigma_plus x $tmp, so $tmp is now SIGMA^+ x $svd_u x $y
    $tmp /= $svd_w->setvaltobad(0.0);
    $tmp->inplace->setbadtoval(0.0);

    #EXPLANATION:   $svd_v x SIGMA^+ x $svd_u x $y
    return sumover( $svd_v * $tmp );

} # sub: _svd()

#_mkbasis returns an ndarray in which the k(=j*n+i)_th column is v**j * u**i
#k=0 j=0 i=0
#k=1 j=0 i=1
#k=2 j=0 i=2
#k=3 j=1 i=0
# ...

#each row corresponds to a control point
#and the rows for each of the control points look like this, e.g.:
# (1) (u) (u**2) (v) (vu) (v(u**2)) (v**2) ((v**2)u) ((v**2)(u**2))
# and so on for the next control point.

sub _mkbasis ($$$$) {
    my $fit    = shift; # dims n n
    my $npts   = shift; # scalar
    my $u      = shift; # dims npts
    my $v      = shift; # dims npts
    my $ncoeff = sum( $fit );
    my $fit_coords = $fit->whichND; # dims uv ncoeff
    cat($u,$v)            # npts uv
      ->transpose         # uv npts
      ->dummy(1,$ncoeff)  # uv ncoeff npts
      ->ipow($fit_coords) # same
      ->prodover          # ncoeff npts
      ;
} # sub: _mkbasis()

sub PDL::fitwarp2d {
    croak "Usage: (\$px,\$py) = fitwarp2d(x(m);y(m);u(m);v(m);\$nf; { options })"
	if $#_ < 4 or ( $#_ >= 5 and ref($_[5]) ne "HASH" );

    my $x  = shift;
    my $y  = shift;
    my $u  = shift;
    my $v  = shift;
    my $nf = shift;

    my $opts = PDL::Options->new( { FIT => ones(byte,$nf,$nf) } ); #, THRESH => 1.0e-5 } );
    $opts->options( $_[0] ) if $#_ > -1;
    my $oref = $opts->current();

    # safety checks
    my $npts = $x->nelem;
    croak "fitwarp2d: x, y, u, and v must be the same size (and 1D)"
	unless $npts == $y->nelem and $npts == $u->nelem and $npts == $v->nelem
	    and $x->getndims == 1 and $y->getndims == 1 and $u->getndims == 1 and $v->getndims == 1;

#    my $svd_thresh = $$oref{THRESH};
#    croak "fitwarp2d: THRESH option must be >= 0."
#	if $svd_thresh < 0;

    my $fit = $$oref{FIT};
    my $fit_ndim = $fit->getndims();
    croak "fitwarp2d: FIT option must be sent a (\$nf,\$nf[,2]) element ndarray"
	unless UNIVERSAL::isa($fit,"PDL") and
	    ($fit_ndim == 2 or ($fit_ndim == 3 and $fit->getdim(2) == 2)) and
	    $fit->getdim(0) == $nf and $fit->getdim(1) == $nf;

    # how many coeffs to fit (first we ensure $fit is either 0 or 1)
    $fit = convert( $fit != 0, byte );

    my ( $fitx, $fity, $ncoeffx, $ncoeffy, $ncoeff );
    if ( $fit_ndim == 2 ) {
	$fitx = $fit;
	$fity = $fit;
	$ncoeff = $ncoeffx = $ncoeffy = sum( $fit );
    } else {
	$fitx = $fit->slice(",,(0)");
	$fity = $fit->slice(",,(1)");
	$ncoeffx = sum($fitx);
	$ncoeffy = sum($fity);
	$ncoeff = $ncoeffx > $ncoeffy ? $ncoeffx : $ncoeffy;
    }

    croak "fitwarp2d: number of points ($npts) must be >= \$ncoeff ($ncoeff)"
	unless $npts >= $ncoeff;

    # create the basis functions for the SVD fitting
    my ( $basisx, $basisy );

    $basisx = _mkbasis( $fitx, $npts, $u, $v );

    if ( $fit_ndim == 2 ) {
	$basisy = $basisx;
    } else {
	$basisy = _mkbasis( $fity, $npts, $u, $v );
    }

    my $px = _svd( $basisx, $x ); # $svd_thresh);
    my $py = _svd( $basisy, $y ); # $svd_thresh);

    # convert into $nf x $nf element ndarrays, if necessary
    my $nf2 = $nf * $nf;

    return ( $px->reshape($nf,$nf), $py->reshape($nf,$nf) )
	if $ncoeff == $nf2 and $ncoeffx == $ncoeffy;

    # re-create the matrix
    my $xtmp = zeroes( $nf, $nf );
    my $ytmp = zeroes( $nf, $nf );

    my $kx = 0;
    my $ky = 0;
    foreach my $i ( 0 .. ($nf - 1) ) {
	foreach my $j ( 0 .. ($nf - 1) ) {
	    if ( $fitx->at($i,$j) ) {
		$xtmp->set($i,$j, $px->at($kx) );
		$kx++;
	    }
	    if ( $fity->at($i,$j) ) {
		$ytmp->set($i,$j, $py->at($ky) );
		$ky++;
	    }
	}
    }

    return ( $xtmp, $ytmp )

} # sub: fitwarp2d

*fitwarp2d = \&PDL::fitwarp2d;

sub PDL::applywarp2d {
    # checks
    croak "Usage: (\$x,\$y) = applywarp2d(px(nf,nf);py(nf,nf);u(m);v(m);)"
	if $#_ != 3;

    my $px = shift;
    my $py = shift;
    my $u  = shift;
    my $v  = shift;
    my $npts = $u->nelem;

    # safety check
    croak "applywarp2d: u and v must be the same size (and 1D)"
	unless $npts == $u->nelem and $npts == $v->nelem
	    and $u->getndims == 1 and $v->getndims == 1;

    my $nf  = $px->getdim(0);
    my $nf2 = $nf * $nf;

    # could remove terms with 0 coeff here
    # (would also have to remove them from px/py for
    #  the matrix multiplication below)
    #
    my $mat = _mkbasis( ones(byte,$nf,$nf), $npts, $u, $v );

    my $x = reshape( $mat x $px->flat->transpose(), $npts );
    my $y = reshape( $mat x $py->flat->transpose(), $npts );
    return ( $x, $y );

} # sub: applywarp2d

*applywarp2d = \&PDL::applywarp2d;

' );

## resampling routines taken from v3.6-0 of the Eclipse package
##   http://www.eso.org/eclipse by Nicolas Devillard
##

pp_addhdr( '#include "resample.h"' . "\n" );

pp_def(
    'warp2d',
    Doc=> <<'EOD',
=for ref

Warp a 2D image given a polynomial describing the I mapping.

=for usage

  $out = warp2d( $img, $px, $py, { options } );

Apply the polynomial transformation encoded in the C<$px> and
C<$py> ndarrays to warp the input image C<$img> into the output
image C<$out>.

The format for the polynomial transformation is described in
the documentation for the L routine.

At each point C, the closest 16 pixel values are combined
with an interpolation kernel to calculate the value at C.
The interpolation is therefore done in the image, rather than
Fourier, domain.
By default, a C kernel is used, but this can be changed
using the C option discussed below
(the choice of kernel depends on the frequency content of the input image).

The routine is based on the C command from
the Eclipse data-reduction package - see http://www.eso.org/eclipse/ - and
for further details on image resampling see
Wolberg, G., "Digital Image Warping", 1990, IEEE Computer
Society Press ISBN 0-8186-8944-7).

Currently the output image is the same size as the input one,
which means data will be lost if the transformation reduces
the pixel scale.  This will (hopefully) be changed soon.

=for example

  $img = rvals(byte,501,501);
  imag $img, { JUSTIFY => 1 };
  #
  # use a not-particularly-obvious transformation:
  #   x = -10 + 0.5 * $u - 0.1 * $v
  #   y = -20 + $v - 0.002 * $u * $v
  #
  $px  = pdl( [ -10, 0.5 ], [ -0.1, 0 ] );
  $py  = pdl( [ -20, 0 ], [ 1, 0.002 ] );
  $wrp = warp2d( $img, $px, $py );
  #
  # see the warped image
  imag $warp, { JUSTIFY => 1 };

The options are:

=for options

  KERNEL - default value is tanh
  NOVAL  - default value is 0

C is used to specify which interpolation kernel to use
(to see what these kernels look like, use the
L routine).
The options are:

=over 4

=item tanh

Hyperbolic tangent: the approximation of an ideal box filter by the
product of symmetric tanh functions.

=item sinc

For a correctly sampled signal, the ideal filter in the fourier domain is a rectangle,
which produces a C interpolation kernel in the spatial domain:

  sinc(x) = sin(pi * x) / (pi * x)

However, it is not ideal for the C<4x4> pixel region used here.

=item sinc2

This is the square of the sinc function.

=item lanczos

Although defined differently to the C kernel, the result is very
similar in the spatial domain.  The Lanczos function is defined as

  L(x) = sinc(x) * sinc(x/2)  if abs(x) < 2
       = 0                       otherwise

=item hann

This kernel is derived from the following function:

  H(x) = a + (1-a) * cos(2*pi*x/(N-1))  if abs(x) < 0.5*(N-1)
       = 0                                 otherwise

with C and N currently equal to 2001.

=item hamming

This kernel uses the same C as the Hann filter, but with
C.

=back

C gives the value used to indicate that a pixel in the
output image does not map onto one in the input image.
EOD
    HandleBad => 0,
    Pars => 'img(m,n); ldouble px(np,np); ldouble py(np,np); [o] warp(m,n); ldouble [t] poly(np); ldouble [t] kernel(ns)',
    OtherPars => 'char *kernel_type; double noval; IV nsamples => ns',
    GenericTypes => [ qw(F D E) ],
        PMCode => <<'EOF',
# support routine
{
    my %warp2d = map { ($_,1) } qw( tanh sinc sinc2 lanczos hamming hann );
    # note: convert to lower case
    sub _check_kernel ($$) {
	my $kernel = lc shift;
	my $code   = shift;
	barf "Unknown kernel $kernel sent to $code\n" .
	    "\tmust be one of [" . join(',',keys %warp2d) . "]\n"
		unless exists $warp2d{$kernel};
	return $kernel;
    }
}

sub PDL::warp2d {
    my $opts = PDL::Options->new( { KERNEL => "tanh", NOVAL => 0 } );
    $opts->options( pop(@_) ) if ref($_[$#_]) eq "HASH";
    die "Usage: warp2d( in(m,n), px(np,np); py(np,np); [o] out(m,n), {Options} )"
	if $#_<2 || $#_>3;
    my $img = shift;
    my $px  = shift;
    my $py  = shift;
    my $out = $#_ == -1 ? PDL->null() : shift;
    # safety checks
    my $copt   = $opts->current();
    my $kernel = _check_kernel( $$copt{KERNEL}, "warp2d" );
    &PDL::_warp2d_int( $img, $px, $py, $out, $kernel, $$copt{NOVAL}, _get_kernel_size() );
    return $out;
}
EOF
    Code => <<'EOF',
    PDL_Indx k, lx_3, ly_3, px, py, tabx, taby;
    PDL_Indx da[16], db[16] ;
    PDL_LDouble cur, neighbors[16], rsc[8], sumrs, x, y, *kernel = $P(kernel), *poly = $P(poly);

    /* Generate interpolation kernel */
    if (!generate_interpolation_kernel($COMP(kernel_type), $SIZE(ns), kernel))
      $CROAK("Invalid kernel type '%s'", $COMP(kernel_type));

    /* Compute sizes  */
    lx_3 = $SIZE(m) - 3;
    ly_3 = $SIZE(n) - 3;

    /* Pre compute leaps for 16 closest neighbors positions */
    da[0]  = -1;  db[0]  = -1;
    da[1]  =  0;  db[1]  = -1;
    da[2]  =  1;  db[2]  = -1;
    da[3]  =  2;  db[3]  = -1;

    da[4]  = -1;  db[4]  =  0;
    da[5]  =  0;  db[5]  =  0;
    da[6]  =  1;  db[6]  =  0;
    da[7]  =  2;  db[7]  =  0;

    da[8]  = -1;  db[8]  =  1;
    da[9]  =  0;  db[9]  =  1;
    da[10] =  1;  db[10] =  1;
    da[11] =  2;  db[11] =  1;

    da[12] = -1;  db[12] =  2;
    da[13] =  0;  db[13] =  2;
    da[14] =  1;  db[14] =  2;
    da[15] =  2;  db[15] =  2;

    poly[0] = 1.0;

    /* Loop over the output image */
    broadcastloop %{
	loop(n) %{

	    /* fill in poly array */
	    loop (np=1) %{
		poly[np] = (PDL_LDouble) n * poly[np-1];
	    %}

	    loop(m) %{

		/* Compute the original source for this pixel   */
                x = poly2d_compute( $SIZE(np), $P(px), m, poly );
                y = poly2d_compute( $SIZE(np), $P(py), m, poly );

		/* Which is the closest integer positioned neighbor? */
		px = (int)x ;
		py = (int)y ;

		if ((px < 1) || (px > lx_3) || (py < 1) || (py > ly_3))
		    $warp() = ($GENERIC()) $COMP(noval);
		else {

		    /* Now feed the positions for the closest 16 neighbors  */
		    for (k=0 ; k<16 ; k++) {
			neighbors[k] = (PDL_LDouble) $img( m => px+da[k], n => py+db[k] );
		    }

		    /* Which tabulated value index shall we use?    */
		    tabx = (x - (PDL_LDouble)px) * (PDL_LDouble)(TABSPERPIX) ;
		    taby = (y - (PDL_LDouble)py) * (PDL_LDouble)(TABSPERPIX) ;

		    /* Compute resampling coefficients  */
		    /* rsc[0..3] in x, rsc[4..7] in y   */

		    rsc[0] = kernel[TABSPERPIX + tabx] ;
		    rsc[1] = kernel[tabx] ;
		    rsc[2] = kernel[TABSPERPIX - tabx] ;
		    rsc[3] = kernel[2 * TABSPERPIX - tabx] ;
		    rsc[4] = kernel[TABSPERPIX + taby] ;
		    rsc[5] = kernel[taby] ;
		    rsc[6] = kernel[TABSPERPIX - taby] ;
		    rsc[7] = kernel[2 * TABSPERPIX - taby] ;

		    sumrs = (rsc[0]+rsc[1]+rsc[2]+rsc[3]) *
			(rsc[4]+rsc[5]+rsc[6]+rsc[7]) ;

		    /* Compute interpolated pixel now   */
		    cur =   rsc[4] * (  rsc[0]*neighbors[0] +
					rsc[1]*neighbors[1] +
					rsc[2]*neighbors[2] +
					rsc[3]*neighbors[3] ) +
			    rsc[5] * (  rsc[0]*neighbors[4] +
					rsc[1]*neighbors[5] +
					rsc[2]*neighbors[6] +
					rsc[3]*neighbors[7] ) +
			    rsc[6] * (  rsc[0]*neighbors[8] +
					rsc[1]*neighbors[9] +
					rsc[2]*neighbors[10] +
					rsc[3]*neighbors[11] ) +
			    rsc[7] * (  rsc[0]*neighbors[12] +
					rsc[1]*neighbors[13] +
					rsc[2]*neighbors[14] +
					rsc[3]*neighbors[15] ) ;

		    /* Copy the value to the output image */
		    $warp() = ($GENERIC()) (cur/sumrs);

		} /* if: edge or interior */

	    %}   /* loop(m) */
	 %}   /* loop(n) */
    %}        /* broadcastloop */
EOF
    ); # pp_def: warp2d

pp_addxs( '

int
_get_kernel_size()
  PROTOTYPE:
  CODE:
    RETVAL = KERNEL_SAMPLES;
  OUTPUT:
    RETVAL

');

pp_def( 'warp2d_kernel',
    Doc => <<'EOF',
=for ref

Return the specified kernel, as used by L

=for usage

  ( $x, $k ) = warp2d_kernel( $name )

The valid values for C<$name> are the same as the C option
of L.

=for example

  line warp2d_kernel( "hamming" );
EOF
    HandleBad => 0,
    PMCode => '
sub PDL::warp2d_kernel ($) {
    my $kernel = _check_kernel( shift, "warp2d_kernel" );
    &PDL::_warp2d_kernel_int( my $x=PDL->null, my $k=PDL->null, $kernel, _get_kernel_size() );
    return ( $x, $k );
}
',
     Pars => '[o] x(n); [o] k(n); ldouble [t] kernel(n)',
     OtherPars => 'char *name; PDL_Indx nsize => n',
     GenericTypes => [ qw(F D E) ],
     Code => <<'EOF',
    PDL_LDouble     *kernel = $P(kernel);
    if (!generate_interpolation_kernel($COMP(name), $SIZE(n), kernel))
      $CROAK("Invalid kernel type '%s'", $COMP(name));
    /* fill in ndarrays */
    PDL_LDouble xx = 0.0;
    broadcastloop %{
	loop (n) %{
	    $x() = xx;
	    $k() = kernel[n];
	    xx += 1.0 / (double) TABSPERPIX;
	%}
    %}
EOF
);

pp_done();