"Fossies" - the Fresh Open Source Software Archive

Member "PDL-2.080/GENERATED/PDL/Image2D.pm" (28 May 2022, 37883 Bytes) of package /linux/misc/PDL-2.080.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Perl source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "Image2D.pm" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 2.079_vs_2.080.

    1 #
    2 # GENERATED WITH PDL::PP! Don't modify!
    3 #
    4 package PDL::Image2D;
    5 
    6 our @EXPORT_OK = qw( conv2d med2d med2df box2d patch2d patchbad2d max2d_ind centroid2d crop cc8compt cc4compt ccNcompt polyfill pnpoly polyfillv rotnewsz rot2d bilin2d rescale2d fitwarp2d applywarp2d warp2d warp2d_kernel warp2d_kernel );
    7 our %EXPORT_TAGS = (Func=>\@EXPORT_OK);
    8 
    9 use PDL::Core;
   10 use PDL::Exporter;
   11 use DynaLoader;
   12 
   13 
   14    
   15    our @ISA = ( 'PDL::Exporter','DynaLoader' );
   16    push @PDL::Core::PP, __PACKAGE__;
   17    bootstrap PDL::Image2D ;
   18 
   19 
   20 
   21 
   22 
   23 
   24 #line 4 "image2d.pd"
   25 
   26 use strict;
   27 use warnings;
   28 
   29 =head1 NAME
   30 
   31 PDL::Image2D - Miscellaneous 2D image processing functions
   32 
   33 =head1 DESCRIPTION
   34 
   35 Miscellaneous 2D image processing functions - for want
   36 of anywhere else to put them.
   37 
   38 =head1 SYNOPSIS
   39 
   40  use PDL::Image2D;
   41 
   42 =cut
   43 
   44 use PDL;  # ensure qsort routine available
   45 use PDL::Math;
   46 use Carp;
   47 #line 48 "Image2D.pm"
   48 
   49 
   50 
   51 
   52 
   53 
   54 =head1 FUNCTIONS
   55 
   56 =cut
   57 
   58 
   59 
   60 
   61 #line 950 "../../blib/lib/PDL/PP.pm"
   62 #line 63 "Image2D.pm"
   63 
   64 
   65 
   66 #line 950 "../../blib/lib/PDL/PP.pm"
   67 #line 68 "Image2D.pm"
   68 
   69 
   70 
   71 #line 950 "../../blib/lib/PDL/PP.pm"
   72 #line 73 "Image2D.pm"
   73 
   74 
   75 
   76 #line 948 "../../blib/lib/PDL/PP.pm"
   77 
   78 
   79 
   80 =head2 conv2d
   81 
   82 =for sig
   83 
   84   Signature: (a(m,n); kern(p,q); [o]b(m,n); int opt)
   85 
   86 =for ref
   87 
   88 2D convolution of an array with a kernel (smoothing)
   89 
   90 For large kernels, using a FFT routine,
   91 such as L<fftconvolve()|PDL::FFT/fftconvolve()> in C<PDL::FFT>,
   92 will be quicker.
   93 
   94 =for usage
   95 
   96  $new = conv2d $old, $kernel, {OPTIONS}
   97 
   98 =for example
   99 
  100  $smoothed = conv2d $image, ones(3,3), {Boundary => Reflect}
  101 
  102 =for options
  103 
  104  Boundary - controls what values are assumed for the image when kernel
  105             crosses its edge:
  106         => Default   - periodic boundary conditions
  107                            (i.e. wrap around axis)
  108         => Reflect   - reflect at boundary
  109         => Truncate  - truncate at boundary
  110         => Replicate - repeat boundary pixel values
  111 
  112 
  113 
  114 =for bad
  115 
  116 Unlike the FFT routines, conv2d is able to process bad values.
  117 
  118 =cut
  119 #line 120 "Image2D.pm"
  120 
  121 
  122 
  123 #line 949 "../../blib/lib/PDL/PP.pm"
  124 
  125 
  126 sub PDL::conv2d {
  127    my $opt; $opt = pop @_ if ref($_[$#_]) eq 'HASH';
  128    die 'Usage: conv2d( a(m,n), kern(p,q), [o]b(m,n), {Options} )'
  129       if $#_<1 || $#_>2;
  130    my($x,$kern) = @_;
  131    my $c = $#_ == 2 ? $_[2] : $x->nullcreate;
  132    PDL::_conv2d_int($x,$kern,$c,
  133     (!(defined $opt && exists $$opt{Boundary}))?0:
  134     (($$opt{Boundary} eq "Reflect") +
  135     2*($$opt{Boundary} eq "Truncate") +
  136     3*($$opt{Boundary} eq "Replicate")));
  137    return $c;
  138 }
  139 #line 140 "Image2D.pm"
  140 
  141 
  142 
  143 #line 950 "../../blib/lib/PDL/PP.pm"
  144 
  145 *conv2d = \&PDL::conv2d;
  146 #line 147 "Image2D.pm"
  147 
  148 
  149 
  150 #line 948 "../../blib/lib/PDL/PP.pm"
  151 
  152 
  153 
  154 =head2 med2d
  155 
  156 =for sig
  157 
  158   Signature: (a(m,n); kern(p,q); [o]b(m,n); double [t]tmp(pq); int opt)
  159 
  160 =for ref
  161 
  162 2D median-convolution of an array with a kernel (smoothing)
  163 
  164 Note: only points in the kernel E<gt>0 are included in the median, other
  165 points are weighted by the kernel value (medianing lots of zeroes
  166 is rather pointless)
  167 
  168 =for usage
  169 
  170  $new = med2d $old, $kernel, {OPTIONS}
  171 
  172 =for example
  173 
  174  $smoothed = med2d $image, ones(3,3), {Boundary => Reflect}
  175 
  176 =for options
  177 
  178  Boundary - controls what values are assumed for the image when kernel
  179             crosses its edge:
  180         => Default   - periodic boundary conditions (i.e. wrap around axis)
  181         => Reflect   - reflect at boundary
  182         => Truncate  - truncate at boundary
  183         => Replicate - repeat boundary pixel values
  184 
  185 
  186 
  187 =for bad
  188 
  189 Bad values are ignored in the calculation. If all elements within the
  190 kernel are bad, the output is set bad.
  191 
  192 =cut
  193 #line 194 "Image2D.pm"
  194 
  195 
  196 
  197 #line 949 "../../blib/lib/PDL/PP.pm"
  198 
  199 sub PDL::med2d {
  200    my $opt; $opt = pop @_ if ref($_[$#_]) eq 'HASH';
  201    die 'Usage: med2d( a(m,n), kern(p,q), [o]b(m,n), {Options} )'
  202       if $#_<1 || $#_>2;
  203    my($x,$kern) = @_;
  204    croak "med2d: kernel must contain some positive elements.\n"
  205        if all( $kern <= 0 );
  206    my $c = $#_ == 2 ? $_[2] : $x->nullcreate;
  207    PDL::_med2d_int($x,$kern,$c,
  208     (!(defined $opt && exists $$opt{Boundary}))?0:
  209     (($$opt{Boundary} eq "Reflect") +
  210     2*($$opt{Boundary} eq "Truncate") +
  211     3*($$opt{Boundary} eq "Replicate")));
  212    return $c;
  213 }
  214 #line 215 "Image2D.pm"
  215 
  216 
  217 
  218 #line 950 "../../blib/lib/PDL/PP.pm"
  219 
  220 *med2d = \&PDL::med2d;
  221 #line 222 "Image2D.pm"
  222 
  223 
  224 
  225 #line 948 "../../blib/lib/PDL/PP.pm"
  226 
  227 
  228 
  229 =head2 med2df
  230 
  231 =for sig
  232 
  233   Signature: (a(m,n); [o]b(m,n); int __p_size; int __q_size; int opt)
  234 
  235 
  236 =for ref
  237 
  238 2D median-convolution of an array in a pxq window (smoothing)
  239 
  240 Note: this routine does the median over all points in a rectangular
  241       window and is not quite as flexible as C<med2d> in this regard
  242       but slightly faster instead
  243 
  244 =for usage
  245 
  246  $new = med2df $old, $xwidth, $ywidth, {OPTIONS}
  247 
  248 =for example
  249 
  250  $smoothed = med2df $image, 3, 3, {Boundary => Reflect}
  251 
  252 =for options
  253 
  254  Boundary - controls what values are assumed for the image when kernel
  255             crosses its edge:
  256         => Default   - periodic boundary conditions (i.e. wrap around axis)
  257         => Reflect   - reflect at boundary
  258         => Truncate  - truncate at boundary
  259         => Replicate - repeat boundary pixel values
  260 
  261 
  262 
  263 =for bad
  264 
  265 med2df does not process bad values.
  266 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  267 
  268 
  269 =cut
  270 #line 271 "Image2D.pm"
  271 
  272 
  273 
  274 #line 949 "../../blib/lib/PDL/PP.pm"
  275 
  276 sub PDL::med2df {
  277    my $opt; $opt = pop @_ if ref($_[$#_]) eq 'HASH';
  278    die 'Usage: med2df( a(m,n), [o]b(m,n), p, q, {Options} )'
  279       if $#_<2 || $#_>3;
  280    my($x,$p,$q) = @_;
  281    croak "med2df: kernel must contain some positive elements.\n"
  282        if $p == 0 && $q == 0;
  283    my $c = $#_ == 3 ? $_[3] : $x->nullcreate;
  284    &PDL::_med2df_int($x,$c,$p,$q,
  285     (!(defined $opt && exists $$opt{Boundary}))?0:
  286     (($$opt{Boundary} eq "Reflect") +
  287     2*($$opt{Boundary} eq "Truncate") +
  288     3*($$opt{Boundary} eq "Replicate")));
  289    return $c;
  290 }
  291 #line 292 "Image2D.pm"
  292 
  293 
  294 
  295 #line 950 "../../blib/lib/PDL/PP.pm"
  296 
  297 *med2df = \&PDL::med2df;
  298 #line 299 "Image2D.pm"
  299 
  300 
  301 
  302 #line 948 "../../blib/lib/PDL/PP.pm"
  303 
  304 
  305 
  306 =head2 box2d
  307 
  308 =for sig
  309 
  310   Signature: (a(n,m); [o] b(n,m); int wx; int wy; int edgezero)
  311 
  312 
  313 =for ref
  314 
  315 fast 2D boxcar average
  316 
  317 =for example
  318 
  319   $smoothim = $im->box2d($wx,$wy,$edgezero=1);
  320 
  321 The edgezero argument controls if edge is set to zero (edgezero=1)
  322 or just keeps the original (unfiltered) values.
  323 
  324 C<box2d> should be updated to support similar edge options
  325 as C<conv2d> and C<med2d> etc.
  326 
  327 Boxcar averaging is a pretty crude way of filtering. For serious stuff
  328 better filters are around (e.g., use L</conv2d> with the appropriate
  329 kernel). On the other hand it is fast and computational cost grows only
  330 approximately linearly with window size.
  331 
  332 
  333 
  334 =for bad
  335 
  336 box2d does not process bad values.
  337 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  338 
  339 
  340 =cut
  341 #line 342 "Image2D.pm"
  342 
  343 
  344 
  345 #line 950 "../../blib/lib/PDL/PP.pm"
  346 
  347 *box2d = \&PDL::box2d;
  348 #line 349 "Image2D.pm"
  349 
  350 
  351 
  352 #line 948 "../../blib/lib/PDL/PP.pm"
  353 
  354 
  355 
  356 =head2 patch2d
  357 
  358 =for sig
  359 
  360   Signature: (a(m,n); int bad(m,n); [o]b(m,n))
  361 
  362 
  363 =for ref
  364 
  365 patch bad pixels out of 2D images using a mask
  366 
  367 =for usage
  368 
  369  $patched = patch2d $data, $bad;
  370 
  371 C<$bad> is a 2D mask array where 1=bad pixel 0=good pixel.
  372 Pixels are replaced by the average of their non-bad neighbours;
  373 if all neighbours are bad, the original data value is
  374 copied across.
  375 
  376 
  377 
  378 =for bad
  379 
  380 This routine does not handle bad values - use L</patchbad2d> instead
  381 
  382 =cut
  383 #line 384 "Image2D.pm"
  384 
  385 
  386 
  387 #line 950 "../../blib/lib/PDL/PP.pm"
  388 
  389 *patch2d = \&PDL::patch2d;
  390 #line 391 "Image2D.pm"
  391 
  392 
  393 
  394 #line 948 "../../blib/lib/PDL/PP.pm"
  395 
  396 
  397 
  398 =head2 patchbad2d
  399 
  400 =for sig
  401 
  402   Signature: (a(m,n); [o]b(m,n))
  403 
  404 =for ref
  405 
  406 patch bad pixels out of 2D images containing bad values
  407 
  408 =for usage
  409 
  410  $patched = patchbad2d $data;
  411 
  412 Pixels are replaced by the average of their non-bad neighbours;
  413 if all neighbours are bad, the output is set bad.
  414 If the input ndarray contains I<no> bad values, then a straight copy
  415 is performed (see L</patch2d>).
  416 
  417 
  418 =for bad
  419 
  420 patchbad2d handles bad values. The output ndarray I<may> contain
  421 bad values, depending on the pattern of bad values in the input ndarray.
  422 
  423 =cut
  424 #line 425 "Image2D.pm"
  425 
  426 
  427 
  428 #line 950 "../../blib/lib/PDL/PP.pm"
  429 
  430 *patchbad2d = \&PDL::patchbad2d;
  431 #line 432 "Image2D.pm"
  432 
  433 
  434 
  435 #line 948 "../../blib/lib/PDL/PP.pm"
  436 
  437 
  438 
  439 =head2 max2d_ind
  440 
  441 =for sig
  442 
  443   Signature: (a(m,n); [o]val(); int [o]x(); int[o]y())
  444 
  445 
  446 =for ref
  447 
  448 Return value/position of maximum value in 2D image
  449 
  450 Contributed by Tim Jeness
  451 
  452 
  453 
  454 =for bad
  455 
  456 Bad values are excluded from the search. If all pixels
  457 are bad then the output is set bad.
  458 
  459 
  460 
  461 =cut
  462 #line 463 "Image2D.pm"
  463 
  464 
  465 
  466 #line 950 "../../blib/lib/PDL/PP.pm"
  467 
  468 *max2d_ind = \&PDL::max2d_ind;
  469 #line 470 "Image2D.pm"
  470 
  471 
  472 
  473 #line 948 "../../blib/lib/PDL/PP.pm"
  474 
  475 
  476 
  477 =head2 centroid2d
  478 
  479 =for sig
  480 
  481   Signature: (im(m,n); x(); y(); box(); [o]xcen(); [o]ycen())
  482 
  483 =for ref
  484 
  485 Refine a list of object positions in 2D image by centroiding in a box
  486 
  487 C<$box> is the full-width of the box, i.e. the window
  488 is C<+/- $box/2>.
  489 
  490 
  491 
  492 =for bad
  493 
  494 Bad pixels are excluded from the centroid calculation. If all elements are
  495 bad (or the pixel sum is 0 - but why would you be centroiding
  496 something with negatives in...) then the output values are set bad.
  497 
  498 
  499 =cut
  500 #line 501 "Image2D.pm"
  501 
  502 
  503 
  504 #line 950 "../../blib/lib/PDL/PP.pm"
  505 
  506 *centroid2d = \&PDL::centroid2d;
  507 #line 508 "Image2D.pm"
  508 
  509 
  510 
  511 #line 797 "image2d.pd"
  512 
  513 =head2 crop
  514 
  515 =for ref
  516 
  517 Return bounding box of given mask in an C<indx> ndarray, so it can broadcast.
  518 Use other operations (such as L<PDL::Bad/isgood>, or
  519 L<PDL::Primitive/eqvec> with a colour vector) to create a mask suitable
  520 for your application.
  521 
  522 =for example
  523 
  524   $x1x2y1y2 = crop($image);
  525 
  526 =cut
  527 
  528 *crop = \&PDL::crop;
  529 sub PDL::crop {
  530   my ($mask) = @_;
  531   $mask->xchg(0,1)->orover->_which_int(my $out = null, null);
  532   $out->badflag(1); $out->badvalue(-1);
  533   my ($x1, $x2) = $out->minmaximum;
  534   $mask->orover->_which_int($out = null, null);
  535   $out->badflag(1); $out->badvalue(-1);
  536   my ($y1, $y2) = $out->minmaximum;
  537   $x1->cat($x2, $y1, $y2)->mv(-1,0);
  538 }
  539 #line 540 "Image2D.pm"
  540 
  541 
  542 
  543 #line 827 "image2d.pd"
  544 
  545 
  546 =head2 cc8compt
  547 
  548 =for ref
  549 
  550 Connected 8-component labeling of a binary image.
  551 
  552 Connected 8-component labeling of 0,1 image - i.e. find separate
  553 segmented objects and fill object pixels with object number.
  554 8-component labeling includes all neighboring pixels.
  555 This is just a front-end to ccNcompt.  See also L</cc4compt>.
  556 
  557 =for example
  558 
  559  $segmented = cc8compt( $image > $threshold );
  560 
  561 =head2 cc4compt
  562 
  563 =for ref
  564 
  565 Connected 4-component labeling of a binary image.
  566 
  567 Connected 4-component labeling of 0,1 image - i.e. find separate
  568 segmented objects and fill object pixels with object number.
  569 4-component labling does not include the diagonal neighbors.
  570 This is just a front-end to ccNcompt.  See also L</cc8compt>.
  571 
  572 =for example
  573 
  574  $segmented = cc4compt( $image > $threshold );
  575 
  576 =cut
  577 
  578 sub PDL::cc8compt{
  579 return ccNcompt(shift,8);
  580 }
  581 *cc8compt = \&PDL::cc8compt;
  582 
  583 sub PDL::cc4compt{
  584 return ccNcompt(shift,4);
  585 }
  586 *cc4compt = \&PDL::cc4compt;
  587 #line 588 "Image2D.pm"
  588 
  589 
  590 
  591 #line 948 "../../blib/lib/PDL/PP.pm"
  592 
  593 
  594 
  595 =head2 ccNcompt
  596 
  597 =for sig
  598 
  599   Signature: (a(m,n); int+ [o]b(m,n); int con)
  600 
  601 
  602 
  603 =for ref
  604 
  605 Connected component labeling of a binary image.
  606 
  607 Connected component labeling of 0,1 image - i.e. find separate
  608 segmented objects and fill object pixels with object number.
  609 See also L</cc4compt> and L</cc8compt>.
  610 
  611 The connectivity parameter must be 4 or 8.
  612 
  613 =for example
  614 
  615  $segmented = ccNcompt( $image > $threshold, 4);
  616 
  617  $segmented2 = ccNcompt( $image > $threshold, 8);
  618 
  619 where the second parameter specifies the connectivity (4 or 8) of the labeling.
  620 
  621 
  622 
  623 =for bad
  624 
  625 ccNcompt ignores the bad-value flag of the input ndarrays.
  626 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  627 
  628 
  629 =cut
  630 #line 631 "Image2D.pm"
  631 
  632 
  633 
  634 #line 950 "../../blib/lib/PDL/PP.pm"
  635 
  636 *ccNcompt = \&PDL::ccNcompt;
  637 #line 638 "Image2D.pm"
  638 
  639 
  640 
  641 #line 996 "image2d.pd"
  642 
  643 =head2 polyfill
  644 
  645 =for ref
  646 
  647 fill the area of the given polygon with the given colour.
  648 
  649 This function works inplace, i.e. modifies C<im>.
  650 
  651 =for usage
  652 
  653   polyfill($im,$ps,$colour,[\%options]);
  654 
  655 The default method of determining which points lie inside of the polygon used
  656 is not as strict as the method used in L</pnpoly>. Often, it includes vertices
  657 and edge points. Set the C<Method> option to change this behaviour.
  658 
  659 =for option
  660 
  661 Method   -  Set the method used to determine which points lie in the polygon.
  662             => Default - internal PDL algorithm
  663             => pnpoly  - use the L</pnpoly> algorithm
  664 
  665 =for example
  666 
  667   # Make a convex 3x3 square of 1s in an image using the pnpoly algorithm
  668   $ps = pdl([3,3],[3,6],[6,6],[6,3]);
  669   polyfill($im,$ps,1,{'Method' =>'pnpoly'});
  670 
  671 =cut
  672 sub PDL::polyfill {
  673     my $opt;
  674     $opt = pop @_ if ref($_[-1]) eq 'HASH';
  675     barf('Usage: polyfill($im,$ps,$colour,[\%options])') unless @_==3;
  676     my ($im,$ps,$colour) = @_;
  677 
  678     if($opt) {
  679         use PDL::Options qw();
  680         my $parsed = PDL::Options->new({'Method' => undef});
  681         $parsed->options($opt);
  682         if( $parsed->current->{'Method'} eq 'pnpoly' ) {
  683             PDL::pnpolyfill_pp($im,$ps,$colour);
  684         }
  685     }
  686     else
  687     {
  688         PDL::polyfill_pp($im,$ps,$colour);
  689     }
  690     return $im;
  691 
  692 }
  693 
  694 *polyfill = \&PDL::polyfill;
  695 #line 696 "Image2D.pm"
  696 
  697 
  698 
  699 #line 1053 "image2d.pd"
  700 
  701 
  702 =head2 pnpoly
  703 
  704 =for ref
  705 
  706 'points in a polygon' selection from a 2-D ndarray
  707 
  708 =for usage
  709 
  710   $mask = $img->pnpoly($ps);
  711 
  712   # Old style, do not use
  713   $mask = pnpoly($x, $y, $px, $py);
  714 
  715 For a closed polygon determined by the sequence of points in {$px,$py}
  716 the output of pnpoly is a mask corresponding to whether or not each
  717 coordinate (x,y) in the set of test points, {$x,$y}, is in the interior
  718 of the polygon.  This is the 'points in a polygon' algorithm from
  719 L<http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html>
  720 and vectorized for PDL by Karl Glazebrook.
  721 
  722 =for example
  723 
  724   # define a 3-sided polygon (a triangle)
  725   $ps = pdl([3, 3], [20, 20], [34, 3]);
  726 
  727   # $tri is 0 everywhere except for points in polygon interior
  728   $tri = $img->pnpoly($ps);
  729 
  730   With the second form, the x and y coordinates must also be specified.
  731   B< I<THIS IS MAINTAINED FOR BACKWARD COMPATIBILITY ONLY> >.
  732 
  733   $px = pdl( 3, 20, 34 );
  734   $py = pdl( 3, 20,  3 );
  735   $x = $img->xvals;      # get x pixel coords
  736   $y = $img->yvals;      # get y pixel coords
  737 
  738   # $tri is 0 everywhere except for points in polygon interior
  739   $tri = pnpoly($x,$y,$px,$py);
  740 
  741 =cut
  742 
  743 # From: http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
  744 #
  745 # Fixes needed to pnpoly code:
  746 #
  747 # Use topdl() to ensure ndarray args
  748 #
  749 # Add POD docs for usage
  750 #
  751 # Calculate first term in & expression to use to fix divide-by-zero when
  752 #   the test point is in line with a vertical edge of the polygon.
  753 #   By adding the value of $mask we prevent a divide-by-zero since the &
  754 #   operator does not "short circuit".
  755 
  756 sub PDL::pnpoly {
  757     barf('Usage: $mask = pnpoly($img, $ps);') unless(@_==2 || @_==4);
  758     my ($tx, $ty, $vertx, $verty) = @_;
  759 
  760     # if only two inputs, use the pure PP version
  761     unless( defined $vertx ) {
  762         barf("ps must contain pairwise points.\n") unless $ty->getdim(0) == 2;
  763 
  764         # Input mapping:  $img => $tx, $ps => $ty
  765         return PDL::pnpoly_pp($tx,$ty);
  766     }
  767 
  768     my $testx = PDL::Core::topdl($tx)->dummy(0);
  769     my $testy = PDL::Core::topdl($ty)->dummy(0);
  770     my $vertxj = PDL::Core::topdl($vertx)->rotate(1);
  771     my $vertyj = PDL::Core::topdl($verty)->rotate(1);
  772     my $mask = ( ($verty>$testy) == ($vertyj>$testy) );
  773     my $c = sumover( ! $mask & ( $testx < ($vertxj-$vertx) * ($testy-$verty)
  774                                  / ($vertyj-$verty+$mask) + $vertx) ) % 2;
  775     return $c;
  776 }
  777 
  778 *pnpoly = \&PDL::pnpoly;
  779 #line 780 "Image2D.pm"
  780 
  781 
  782 
  783 #line 1136 "image2d.pd"
  784 
  785 
  786 =head2 polyfillv
  787 
  788 =for ref
  789 
  790 return the (dataflowed) area of an image described by a polygon
  791 
  792 =for usage
  793 
  794   polyfillv($im,$ps,[\%options]);
  795 
  796 The default method of determining which points lie inside of the polygon used
  797 is not as strict as the method used in L</pnpoly>. Often, it includes vertices
  798 and edge points. Set the C<Method> option to change this behaviour.
  799 
  800 =for option
  801 
  802 Method   -  Set the method used to determine which points lie in the polygon.
  803             => Default - internal PDL algorithm
  804             => pnpoly  - use the L</pnpoly> algorithm
  805 
  806 =for example
  807 
  808   # increment intensity in area bounded by $poly using the pnpoly algorithm
  809   $im->polyfillv($poly,{'Method'=>'pnpoly'})++; # legal in perl >= 5.6
  810 
  811   # compute average intensity within area bounded by $poly using the default algorithm
  812   $av = $im->polyfillv($poly)->avg;
  813 
  814 =cut
  815 
  816 sub PDL::polyfillv :lvalue {
  817     my $opt;
  818     $opt = pop @_ if ref($_[-1]) eq 'HASH';
  819     barf('Usage: polyfillv($im,$ps,[\%options])') unless @_==2;
  820 
  821     my ($im,$ps) = @_;
  822     barf("ps must contain pairwise points.\n") unless $ps->getdim(0) == 2;
  823 
  824     if($opt) {
  825         use PDL::Options qw();
  826         my $parsed = PDL::Options->new({'Method' => undef});
  827         $parsed->options($opt);
  828         return $im->where(PDL::pnpoly_pp($im, $ps)) if $parsed->current->{'Method'} eq 'pnpoly';
  829     }
  830 
  831     my $msk = zeroes(long,$im->dims);
  832     PDL::polyfill_pp($msk, $ps, 1);
  833     return $im->where($msk);
  834 }
  835 *polyfillv = \&PDL::polyfillv;
  836 #line 837 "Image2D.pm"
  837 
  838 
  839 
  840 #line 948 "../../blib/lib/PDL/PP.pm"
  841 
  842 
  843 
  844 =head2 rot2d
  845 
  846 =for sig
  847 
  848   Signature: (im(m,n); float angle(); bg(); int aa(); [o] om(p,q))
  849 
  850 
  851 =for ref
  852 
  853 rotate an image by given C<angle>
  854 
  855 =for example
  856 
  857   # rotate by 10.5 degrees with antialiasing, set missing values to 7
  858   $rot = $im->rot2d(10.5,7,1);
  859 
  860 This function rotates an image through an C<angle> between -90 and + 90
  861 degrees. Uses/doesn't use antialiasing depending on the C<aa> flag.
  862 Pixels outside the rotated image are set to C<bg>.
  863 
  864 Code modified from pnmrotate (Copyright Jef Poskanzer) with an algorithm based
  865 on "A Fast Algorithm for General  Raster  Rotation"  by  Alan Paeth,
  866 Graphics Interface '86, pp. 77-81.
  867 
  868 Use the C<rotnewsz> function to find out about the dimension of the
  869 newly created image
  870 
  871   ($newcols,$newrows) = rotnewsz $oldn, $oldm, $angle;
  872 
  873 L<PDL::Transform> offers a more general interface to
  874 distortions, including rotation, with various types of sampling; but
  875 rot2d is faster.
  876 
  877 
  878 
  879 =for bad
  880 
  881 rot2d ignores the bad-value flag of the input ndarrays.
  882 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  883 
  884 
  885 =cut
  886 #line 887 "Image2D.pm"
  887 
  888 
  889 
  890 #line 950 "../../blib/lib/PDL/PP.pm"
  891 
  892 *rot2d = \&PDL::rot2d;
  893 #line 894 "Image2D.pm"
  894 
  895 
  896 
  897 #line 948 "../../blib/lib/PDL/PP.pm"
  898 
  899 
  900 
  901 =head2 bilin2d
  902 
  903 =for sig
  904 
  905   Signature: (Int(n,m); O(q,p))
  906 
  907 
  908 =for ref
  909 
  910 Bilinearly maps the first ndarray in the second. The
  911 interpolated values are actually added to the second
  912 ndarray which is supposed to be larger than the first one.
  913 
  914 
  915 
  916 =for bad
  917 
  918 bilin2d ignores the bad-value flag of the input ndarrays.
  919 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  920 
  921 
  922 =cut
  923 #line 924 "Image2D.pm"
  924 
  925 
  926 
  927 #line 950 "../../blib/lib/PDL/PP.pm"
  928 
  929 *bilin2d = \&PDL::bilin2d;
  930 #line 931 "Image2D.pm"
  931 
  932 
  933 
  934 #line 948 "../../blib/lib/PDL/PP.pm"
  935 
  936 
  937 
  938 =head2 rescale2d
  939 
  940 =for sig
  941 
  942   Signature: (Int(m,n); O(p,q))
  943 
  944 
  945 =for ref
  946 
  947 The first ndarray is rescaled to the dimensions of the second
  948 (expanding or meaning values as needed) and then added to it in place.
  949 Nothing useful is returned.
  950 
  951 If you want photometric accuracy or automatic FITS header metadata
  952 tracking, consider using L<PDL::Transform::map|PDL::Transform/map>
  953 instead: it does these things, at some speed penalty compared to
  954 rescale2d.
  955 
  956 
  957 
  958 =for bad
  959 
  960 rescale2d ignores the bad-value flag of the input ndarrays.
  961 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  962 
  963 
  964 =cut
  965 #line 966 "Image2D.pm"
  966 
  967 
  968 
  969 #line 950 "../../blib/lib/PDL/PP.pm"
  970 
  971 *rescale2d = \&PDL::rescale2d;
  972 #line 973 "Image2D.pm"
  973 
  974 
  975 
  976 #line 1450 "image2d.pd"
  977 
  978 
  979 
  980 =head2 fitwarp2d
  981 
  982 =for ref
  983 
  984 Find the best-fit 2D polynomial to describe
  985 a coordinate transformation.
  986 
  987 =for usage
  988 
  989   ( $px, $py ) = fitwarp2d( $x, $y, $u, $v, $nf, { options } )
  990 
  991 Given a set of points in the output plane (C<$u,$v>), find
  992 the best-fit (using singular-value decomposition) 2D polynomial
  993 to describe the mapping back to the image plane (C<$x,$y>).
  994 The order of the fit is controlled by the C<$nf> parameter
  995 (the maximum power of the polynomial is C<$nf - 1>), and you
  996 can restrict the terms to fit using the C<FIT> option.
  997 
  998 C<$px> and C<$py> are C<np> by C<np> element ndarrays which describe
  999 a polynomial mapping (of order C<np-1>)
 1000 from the I<output> C<(u,v)> image to the I<input> C<(x,y)> image:
 1001 
 1002   x = sum(j=0,np-1) sum(i=0,np-1) px(i,j) * u^i * v^j
 1003   y = sum(j=0,np-1) sum(i=0,np-1) py(i,j) * u^i * v^j
 1004 
 1005 The transformation is returned for the reverse direction (ie
 1006 output to input image) since that is what is required by the
 1007 L<warp2d()|/warp2d> routine.  The L<applywarp2d()|/applywarp2d>
 1008 routine can be used to convert a set of C<$u,$v> points given
 1009 C<$px> and C<$py>.
 1010 
 1011 Options:
 1012 
 1013 =for options
 1014 
 1015   FIT     - which terms to fit? default ones(byte,$nf,$nf)
 1016 
 1017 =begin comment
 1018 
 1019 old option, caused trouble
 1020   THRESH  - in svd, remove terms smaller than THRESH * max value
 1021             default is 1.0e-5
 1022 
 1023 =end comment
 1024 
 1025 =over 4
 1026 
 1027 =item FIT
 1028 
 1029 C<FIT> allows you to restrict which terms of the polynomial to fit:
 1030 only those terms for which the FIT ndarray evaluates to true will be
 1031 evaluated.  If a 2D ndarray is sent in, then it is
 1032 used for the x and y polynomials; otherwise
 1033 C<< $fit->slice(":,:,(0)") >> will be used for C<$px> and
 1034 C<< $fit->slice(":,:,(1)") >> will be used for C<$py>.
 1035 
 1036 =begin comment
 1037 
 1038 =item THRESH
 1039 
 1040 Remove all singular values whose value is less than C<THRESH>
 1041 times the largest singular value.
 1042 
 1043 =end comment
 1044 
 1045 =back
 1046 
 1047 The number of points must be at least equal to the number of
 1048 terms to fit (C<$nf*$nf> points for the default value of C<FIT>).
 1049 
 1050 =for example
 1051 
 1052   # points in original image
 1053   $x = pdl( 0,   0, 100, 100 );
 1054   $y = pdl( 0, 100, 100,   0 );
 1055   # get warped to these positions
 1056   $u = pdl( 10, 10, 90, 90 );
 1057   $v = pdl( 10, 90, 90, 10 );
 1058   #
 1059   # shift of origin + scale x/y axis only
 1060   $fit = byte( [ [1,1], [0,0] ], [ [1,0], [1,0] ] );
 1061   ( $px, $py ) = fitwarp2d( $x, $y, $u, $v, 2, { FIT => $fit } );
 1062   print "px = ${px}py = $py";
 1063   px =
 1064   [
 1065    [-12.5  1.25]
 1066    [    0     0]
 1067   ]
 1068   py =
 1069   [
 1070    [-12.5     0]
 1071    [ 1.25     0]
 1072   ]
 1073   #
 1074   # Compared to allowing all 4 terms
 1075   ( $px, $py ) = fitwarp2d( $x, $y, $u, $v, 2 );
 1076   print "px = ${px}py = $py";
 1077   px =
 1078   [
 1079    [         -12.5           1.25]
 1080    [  1.110223e-16 -1.1275703e-17]
 1081   ]
 1082   py =
 1083   [
 1084    [         -12.5  1.6653345e-16]
 1085    [          1.25 -5.8546917e-18]
 1086   ]
 1087 
 1088   # A higher-degree polynomial should not affect the answer much, but
 1089   # will require more control points
 1090 
 1091   $x = $x->glue(0,pdl(50,12.5, 37.5, 12.5, 37.5));
 1092   $y = $y->glue(0,pdl(50,12.5, 37.5, 37.5, 12.5));
 1093   $u = $u->glue(0,pdl(73,20,40,20,40));
 1094   $v = $v->glue(0,pdl(29,20,40,40,20));
 1095   ( $px3, $py3 ) = fitwarp2d( $x, $y, $u, $v, 3 );
 1096   print "px3 =${px3}py3 =$py3";
 1097   px3 =
 1098   [
 1099    [-6.4981162e+08       71034917     -726498.95]
 1100    [      49902244     -5415096.7      55945.388]
 1101    [    -807778.46      88457.191     -902.51612]
 1102   ]
 1103   py3 =
 1104   [
 1105    [-6.2732159e+08       68576392     -701354.77]
 1106    [      48175125     -5227679.8      54009.114]
 1107    [    -779821.18      85395.681     -871.27997]
 1108   ]
 1109 
 1110   #This illustrates an important point about singular value
 1111   #decompositions that are used in fitwarp2d: like all SVDs, the
 1112   #rotation matrices are not unique, and so the $px and $py returned
 1113   #by fitwarp2d are not guaranteed to be the "simplest" solution.
 1114   #They do still work, though:
 1115 
 1116   ($x3,$y3) = applywarp2d($px3,$py3,$u,$v);
 1117   print approx $x3,$x,1e-4;
 1118   [1 1 1 1 1 1 1 1 1]
 1119   print approx $y3,$y;
 1120   [1 1 1 1 1 1 1 1 1]
 1121 
 1122 =head2 applywarp2d
 1123 
 1124 =for ref
 1125 
 1126 Transform a set of points using a 2-D polynomial mapping
 1127 
 1128 =for usage
 1129 
 1130   ( $x, $y ) = applywarp2d( $px, $py, $u, $v )
 1131 
 1132 Convert a set of points (stored in 1D ndarrays C<$u,$v>)
 1133 to C<$x,$y> using the 2-D polynomial with coefficients stored in C<$px>
 1134 and C<$py>.  See L<fitwarp2d()|/fitwarp2d>
 1135 for more information on the format of C<$px> and C<$py>.
 1136 
 1137 =cut
 1138 
 1139 # use SVD to fit data. Assuming no errors.
 1140 
 1141 =pod
 1142 
 1143 =begin comment
 1144 
 1145 Some explanation of the following three subroutines (_svd, _mkbasis,
 1146 and fitwarp2d): See Wolberg 1990 (full ref elsewhere in this
 1147 documentation), Chapter 3.6 "Polynomial Transformations".  The idea is
 1148 that, given a set of control points in the input and output images
 1149 denoted by coordinates (x,y) and (u,v), we want to create a polynomial
 1150 transformation that relates u to linear combinations of powers of x
 1151 and y, and another that relates v to powers of x and y.
 1152 
 1153 The transformations used here and by Wolberg differ slightly, but the
 1154 basic idea is something like this.  For each u and each v, define a
 1155 transform:
 1156 
 1157 u = (sum over j) (sum over i) a_{ij} x**i * y**j , (eqn 1)
 1158 v = (sum over j) (sum over i) b_{ij} x**i * y**j . (eqn 2)
 1159 
 1160 We can write this in matrix notation.  Given that there are M control
 1161 points, U is a column vector with M rows.  A and B are vectors containing
 1162 the N coefficients (related to the degree of the polynomial fit).  W
 1163 is an MxN matrix of the basis row-vectors (the x**i y**j).
 1164 
 1165 The matrix equations we are trying to solve are
 1166 U = W A , (eqn 3)
 1167 V = W B . (eqn 4)
 1168 
 1169 We need to find the A and B column vectors, those are the coefficients
 1170 of the polynomial terms in W.  W is not square, so it has no inverse.
 1171 But is has a pseudo-inverse W^+ that is NxM.  We are going to use that
 1172 pseudo-inverse to isolate A and B, like so:
 1173 
 1174 W^+ U = W^+ W A = A (eqn 5)
 1175 W^+ V = W^+ W B = B (eqn 6)
 1176 
 1177 We are going to get W^+ by performing a singular value decomposition
 1178 of W (to use some of the variable names below):
 1179 
 1180 W = $svd_u x SIGMA x $svd_v->transpose (eqn 7)
 1181 W^+ = $svd_v x SIGMA^+ x $svd_u->transpose . (eqn 8)
 1182 
 1183 Here SIGMA is a square diagonal matrix that contains the singular
 1184 values of W that are in the variable $svd_w.  SIGMA^+ is the
 1185 pseudo-inverse of SIGMA, which is calculated by replacing the non-zero
 1186 singular values with their reciprocals, and then taking the transpose
 1187 of the matrix (which is a no-op since the matrix is square and
 1188 diagonal).
 1189 
 1190 So the code below does this:
 1191 
 1192 _mkbasis computes the matrix W, given control coordinates u and v and
 1193 the maximum degree of the polynomial (and the terms to use).
 1194 
 1195 _svd takes the svd of W, computes the pseudo-inverse of W, and then
 1196 multiplies that with the U vector (there called $y). The output of
 1197 _svd is the A or B vector in eqns 5 & 6 above. Rarely is the matrix
 1198 multiplication explicit, unfortunately, so I have added EXPLANATIONs
 1199 below.
 1200 
 1201 =end comment
 1202 
 1203 =cut
 1204 
 1205 sub _svd ($$) {
 1206     my $basis  = shift;
 1207     my $y      = shift;
 1208 #    my $thresh = shift;
 1209 
 1210     # if we had errors for these points, would normalise the
 1211     # basis functions, and the output array, by these errors here
 1212 
 1213     # perform the SVD
 1214     my ( $svd_u, $svd_w, $svd_v ) = svd( $basis );
 1215 
 1216     # DAL, 09/2017: the reason for removing ANY singular values, much less
 1217     #the smallest ones, is not clear. I commented the line below since this
 1218     #actually removes the LARGEST values in SIGMA^+.
 1219     # $svd_w *= ( $svd_w >= ($svd_w->max * $thresh ) );
 1220     # The line below would instead remove the SMALLEST values in SIGMA^+, but I can see no reason to include it either.
 1221     # $svd_w *= ( $svd_w <= ($svd_w->min / $thresh ) );
 1222 
 1223     # perform the back substitution
 1224     # EXPLANATION: same thing as $svd_u->transpose x $y->transpose.
 1225     my $tmp = $y x $svd_u;
 1226 
 1227     #EXPLANATION: the division by (the non-zero elements of) $svd_w (the singular values) is
 1228     #equivalent to $sigma_plus x $tmp, so $tmp is now SIGMA^+ x $svd_u x $y
 1229     $tmp /= $svd_w->setvaltobad(0.0);
 1230     $tmp->inplace->setbadtoval(0.0);
 1231 
 1232     #EXPLANATION:   $svd_v x SIGMA^+ x $svd_u x $y
 1233     return sumover( $svd_v * $tmp );
 1234 
 1235 } # sub: _svd()
 1236 
 1237 #_mkbasis returns an ndarray in which the k(=j*n+i)_th column is v**j * u**i
 1238 #k=0 j=0 i=0
 1239 #k=1 j=0 i=1
 1240 #k=2 j=0 i=2
 1241 #k=3 j=1 i=0
 1242 # ...
 1243 
 1244 #each row corresponds to a control point
 1245 #and the rows for each of the control points look like this, e.g.:
 1246 # (1) (u) (u**2) (v) (vu) (v(u**2)) (v**2) ((v**2)u) ((v**2)(u**2))
 1247 # and so on for the next control point.
 1248 
 1249 sub _mkbasis ($$$$) {
 1250     my $fit    = shift;
 1251     my $npts   = shift;
 1252     my $u      = shift;
 1253     my $v      = shift;
 1254 
 1255     my $n      = $fit->getdim(0) - 1;
 1256     my $ncoeff = sum( $fit );
 1257 
 1258     my $basis = zeroes( $u->type, $ncoeff, $npts );
 1259     my $k = 0;
 1260     foreach my $j ( 0 .. $n ) {
 1261     my $tmp_v = $v**$j;
 1262     foreach my $i ( 0 .. $n ) {
 1263         if ( $fit->at($i,$j) ) {
 1264         my $tmp = $basis->slice("($k),:");
 1265         $tmp .= $tmp_v * $u**$i;
 1266         $k++;
 1267         }
 1268     }
 1269     }
 1270     return $basis;
 1271 
 1272 } # sub: _mkbasis()
 1273 
 1274 sub PDL::fitwarp2d {
 1275     croak "Usage: (\$px,\$py) = fitwarp2d(x(m);y(m);u(m);v(m);\$nf; { options })"
 1276     if $#_ < 4 or ( $#_ >= 5 and ref($_[5]) ne "HASH" );
 1277 
 1278     my $x  = shift;
 1279     my $y  = shift;
 1280     my $u  = shift;
 1281     my $v  = shift;
 1282     my $nf = shift;
 1283 
 1284     my $opts = PDL::Options->new( { FIT => ones(byte,$nf,$nf) } ); #, THRESH => 1.0e-5 } );
 1285     $opts->options( $_[0] ) if $#_ > -1;
 1286     my $oref = $opts->current();
 1287 
 1288     # safety checks
 1289     my $npts = $x->nelem;
 1290     croak "fitwarp2d: x, y, u, and v must be the same size (and 1D)"
 1291     unless $npts == $y->nelem and $npts == $u->nelem and $npts == $v->nelem
 1292         and $x->getndims == 1 and $y->getndims == 1 and $u->getndims == 1 and $v->getndims == 1;
 1293 
 1294 #    my $svd_thresh = $$oref{THRESH};
 1295 #    croak "fitwarp2d: THRESH option must be >= 0."
 1296 #   if $svd_thresh < 0;
 1297 
 1298     my $fit = $$oref{FIT};
 1299     my $fit_ndim = $fit->getndims();
 1300     croak "fitwarp2d: FIT option must be sent a (\$nf,\$nf[,2]) element ndarray"
 1301     unless UNIVERSAL::isa($fit,"PDL") and
 1302         ($fit_ndim == 2 or ($fit_ndim == 3 and $fit->getdim(2) == 2)) and
 1303         $fit->getdim(0) == $nf and $fit->getdim(1) == $nf;
 1304 
 1305     # how many coeffs to fit (first we ensure $fit is either 0 or 1)
 1306     $fit = convert( $fit != 0, byte );
 1307 
 1308     my ( $fitx, $fity, $ncoeffx, $ncoeffy, $ncoeff );
 1309     if ( $fit_ndim == 2 ) {
 1310     $fitx = $fit;
 1311     $fity = $fit;
 1312     $ncoeff = $ncoeffx = $ncoeffy = sum( $fit );
 1313     } else {
 1314     $fitx = $fit->slice(",,(0)");
 1315     $fity = $fit->slice(",,(1)");
 1316     $ncoeffx = sum($fitx);
 1317     $ncoeffy = sum($fity);
 1318     $ncoeff = $ncoeffx > $ncoeffy ? $ncoeffx : $ncoeffy;
 1319     }
 1320 
 1321     croak "fitwarp2d: number of points ($npts) must be >= \$ncoeff ($ncoeff)"
 1322     unless $npts >= $ncoeff;
 1323 
 1324     # create the basis functions for the SVD fitting
 1325     my ( $basisx, $basisy );
 1326 
 1327     $basisx = _mkbasis( $fitx, $npts, $u, $v );
 1328 
 1329     if ( $fit_ndim == 2 ) {
 1330     $basisy = $basisx;
 1331     } else {
 1332     $basisy = _mkbasis( $fity, $npts, $u, $v );
 1333     }
 1334 
 1335     my $px = _svd( $basisx, $x ); # $svd_thresh);
 1336     my $py = _svd( $basisy, $y ); # $svd_thresh);
 1337 
 1338     # convert into $nf x $nf element ndarrays, if necessary
 1339     my $nf2 = $nf * $nf;
 1340 
 1341     return ( $px->reshape($nf,$nf), $py->reshape($nf,$nf) )
 1342     if $ncoeff == $nf2 and $ncoeffx == $ncoeffy;
 1343 
 1344     # re-create the matrix
 1345     my $xtmp = zeroes( $nf, $nf );
 1346     my $ytmp = zeroes( $nf, $nf );
 1347 
 1348     my $kx = 0;
 1349     my $ky = 0;
 1350     foreach my $i ( 0 .. ($nf - 1) ) {
 1351     foreach my $j ( 0 .. ($nf - 1) ) {
 1352         if ( $fitx->at($i,$j) ) {
 1353         $xtmp->set($i,$j, $px->at($kx) );
 1354         $kx++;
 1355         }
 1356         if ( $fity->at($i,$j) ) {
 1357         $ytmp->set($i,$j, $py->at($ky) );
 1358         $ky++;
 1359         }
 1360     }
 1361     }
 1362 
 1363     return ( $xtmp, $ytmp )
 1364 
 1365 } # sub: fitwarp2d
 1366 
 1367 *fitwarp2d = \&PDL::fitwarp2d;
 1368 
 1369 sub PDL::applywarp2d {
 1370     # checks
 1371     croak "Usage: (\$x,\$y) = applywarp2d(px(nf,nf);py(nf,nf);u(m);v(m);)"
 1372     if $#_ != 3;
 1373 
 1374     my $px = shift;
 1375     my $py = shift;
 1376     my $u  = shift;
 1377     my $v  = shift;
 1378     my $npts = $u->nelem;
 1379 
 1380     # safety check
 1381     croak "applywarp2d: u and v must be the same size (and 1D)"
 1382     unless $npts == $u->nelem and $npts == $v->nelem
 1383         and $u->getndims == 1 and $v->getndims == 1;
 1384 
 1385     my $nf  = $px->getdim(0);
 1386     my $nf2 = $nf * $nf;
 1387 
 1388     # could remove terms with 0 coeff here
 1389     # (would also have to remove them from px/py for
 1390     #  the matrix multiplication below)
 1391     #
 1392     my $mat = _mkbasis( ones(byte,$nf,$nf), $npts, $u, $v );
 1393 
 1394     my $x = reshape( $mat x $px->clump(-1)->transpose(), $npts );
 1395     my $y = reshape( $mat x $py->clump(-1)->transpose(), $npts );
 1396     return ( $x, $y );
 1397 
 1398 } # sub: applywarp2d
 1399 
 1400 *applywarp2d = \&PDL::applywarp2d;
 1401 #line 1402 "Image2D.pm"
 1402 
 1403 
 1404 
 1405 #line 1885 "image2d.pd"
 1406 
 1407 
 1408 =head2 warp2d
 1409 
 1410 =for sig
 1411 
 1412   Signature: (img(m,n); double px(np,np); double py(np,np); [o] warp(m,n); { options })
 1413 
 1414 =for ref
 1415 
 1416 Warp a 2D image given a polynomial describing the I<reverse> mapping.
 1417 
 1418 =for usage
 1419 
 1420   $out = warp2d( $img, $px, $py, { options } );
 1421 
 1422 Apply the polynomial transformation encoded in the C<$px> and
 1423 C<$py> ndarrays to warp the input image C<$img> into the output
 1424 image C<$out>.
 1425 
 1426 The format for the polynomial transformation is described in
 1427 the documentation for the L<fitwarp2d()|/fitwarp2d> routine.
 1428 
 1429 At each point C<x,y>, the closest 16 pixel values are combined
 1430 with an interpolation kernel to calculate the value at C<u,v>.
 1431 The interpolation is therefore done in the image, rather than
 1432 Fourier, domain.
 1433 By default, a C<tanh> kernel is used, but this can be changed
 1434 using the C<KERNEL> option discussed below
 1435 (the choice of kernel depends on the frequency content of the input image).
 1436 
 1437 The routine is based on the C<warping> command from
 1438 the Eclipse data-reduction package - see http://www.eso.org/eclipse/ - and
 1439 for further details on image resampling see
 1440 Wolberg, G., "Digital Image Warping", 1990, IEEE Computer
 1441 Society Press ISBN 0-8186-8944-7).
 1442 
 1443 Currently the output image is the same size as the input one,
 1444 which means data will be lost if the transformation reduces
 1445 the pixel scale.  This will (hopefully) be changed soon.
 1446 
 1447 =for example
 1448 
 1449   $img = rvals(byte,501,501);
 1450   imag $img, { JUSTIFY => 1 };
 1451   #
 1452   # use a not-particularly-obvious transformation:
 1453   #   x = -10 + 0.5 * $u - 0.1 * $v
 1454   #   y = -20 + $v - 0.002 * $u * $v
 1455   #
 1456   $px  = pdl( [ -10, 0.5 ], [ -0.1, 0 ] );
 1457   $py  = pdl( [ -20, 0 ], [ 1, 0.002 ] );
 1458   $wrp = warp2d( $img, $px, $py );
 1459   #
 1460   # see the warped image
 1461   imag $warp, { JUSTIFY => 1 };
 1462 
 1463 The options are:
 1464 
 1465 =for options
 1466 
 1467   KERNEL - default value is tanh
 1468   NOVAL  - default value is 0
 1469 
 1470 C<KERNEL> is used to specify which interpolation kernel to use
 1471 (to see what these kernels look like, use the
 1472 L<warp2d_kernel()|/warp2d_kernel> routine).
 1473 The options are:
 1474 
 1475 =over 4
 1476 
 1477 =item tanh
 1478 
 1479 Hyperbolic tangent: the approximation of an ideal box filter by the
 1480 product of symmetric tanh functions.
 1481 
 1482 =item sinc
 1483 
 1484 For a correctly sampled signal, the ideal filter in the fourier domain is a rectangle,
 1485 which produces a C<sinc> interpolation kernel in the spatial domain:
 1486 
 1487   sinc(x) = sin(pi * x) / (pi * x)
 1488 
 1489 However, it is not ideal for the C<4x4> pixel region used here.
 1490 
 1491 =item sinc2
 1492 
 1493 This is the square of the sinc function.
 1494 
 1495 =item lanczos
 1496 
 1497 Although defined differently to the C<tanh> kernel, the result is very
 1498 similar in the spatial domain.  The Lanczos function is defined as
 1499 
 1500   L(x) = sinc(x) * sinc(x/2)  if abs(x) < 2
 1501        = 0                       otherwise
 1502 
 1503 =item hann
 1504 
 1505 This kernel is derived from the following function:
 1506 
 1507   H(x) = a + (1-a) * cos(2*pi*x/(N-1))  if abs(x) < 0.5*(N-1)
 1508        = 0                                 otherwise
 1509 
 1510 with C<a = 0.5> and N currently equal to 2001.
 1511 
 1512 =item hamming
 1513 
 1514 This kernel uses the same C<H(x)> as the Hann filter, but with
 1515 C<a = 0.54>.
 1516 
 1517 =back
 1518 
 1519 C<NOVAL> gives the value used to indicate that a pixel in the
 1520 output image does not map onto one in the input image.
 1521 
 1522 =cut
 1523 
 1524 # support routine
 1525 {
 1526     my %warp2d = map { ($_,1) } qw( tanh sinc sinc2 lanczos hamming hann );
 1527 
 1528     # note: convert to lower case
 1529     sub _check_kernel ($$) {
 1530     my $kernel = lc shift;
 1531     my $code   = shift;
 1532     barf "Unknown kernel $kernel sent to $code\n" .
 1533         "\tmust be one of [" . join(',',keys %warp2d) . "]\n"
 1534         unless exists $warp2d{$kernel};
 1535     return $kernel;
 1536     }
 1537 }
 1538 #line 1539 "Image2D.pm"
 1539 
 1540 
 1541 
 1542 #line 949 "../../blib/lib/PDL/PP.pm"
 1543 
 1544 
 1545 
 1546 sub PDL::warp2d {
 1547     my $opts = PDL::Options->new( { KERNEL => "tanh", NOVAL => 0 } );
 1548     $opts->options( pop(@_) ) if ref($_[$#_]) eq "HASH";
 1549 
 1550     die "Usage: warp2d( in(m,n), px(np,np); py(np,np); [o] out(m,n), {Options} )"
 1551     if $#_<2 || $#_>3;
 1552     my $img = shift;
 1553     my $px  = shift;
 1554     my $py  = shift;
 1555     my $out = $#_ == -1 ? PDL->null() : shift;
 1556 
 1557     # safety checks
 1558     my $copt   = $opts->current();
 1559     my $kernel = _check_kernel( $$copt{KERNEL}, "warp2d" );
 1560 
 1561     &PDL::_warp2d_int( $img, $px, $py, $out, $kernel, $$copt{NOVAL} );
 1562     return $out;
 1563 }
 1564 #line 1565 "Image2D.pm"
 1565 
 1566 
 1567 
 1568 #line 950 "../../blib/lib/PDL/PP.pm"
 1569 
 1570 *warp2d = \&PDL::warp2d;
 1571 #line 1572 "Image2D.pm"
 1572 
 1573 
 1574 
 1575 #line 2199 "image2d.pd"
 1576 
 1577 
 1578 
 1579 =head2 warp2d_kernel
 1580 
 1581 =for ref
 1582 
 1583 Return the specified kernel, as used by L</warp2d>
 1584 
 1585 =for usage
 1586 
 1587   ( $x, $k ) = warp2d_kernel( $name )
 1588 
 1589 The valid values for C<$name> are the same as the C<KERNEL> option
 1590 of L<warp2d()|/warp2d>.
 1591 
 1592 =for example
 1593 
 1594   line warp2d_kernel( "hamming" );
 1595 
 1596 =cut
 1597 #line 1598 "Image2D.pm"
 1598 
 1599 
 1600 
 1601 #line 949 "../../blib/lib/PDL/PP.pm"
 1602 
 1603 
 1604 
 1605 sub PDL::warp2d_kernel ($) {
 1606     my $kernel = _check_kernel( shift, "warp2d_kernel" );
 1607 
 1608     my $nelem = _get_kernel_size();
 1609     my $x     = zeroes( $nelem );
 1610     my $k     = zeroes( $nelem );
 1611 
 1612     &PDL::_warp2d_kernel_int( $x, $k, $kernel );
 1613     return ( $x, $k );
 1614 
 1615 #    return _get_kernel( $kernel );
 1616 }
 1617 *warp2d_kernel = \&PDL::warp2d_kernel;
 1618 #line 1619 "Image2D.pm"
 1619 
 1620 
 1621 
 1622 #line 950 "../../blib/lib/PDL/PP.pm"
 1623 
 1624 *warp2d_kernel = \&PDL::warp2d_kernel;
 1625 #line 1626 "Image2D.pm"
 1626 
 1627 
 1628 
 1629 
 1630 
 1631 #line 29 "image2d.pd"
 1632 
 1633 
 1634 =head1 AUTHORS
 1635 
 1636 Copyright (C) Karl Glazebrook 1997 with additions by Robin Williams
 1637 (rjrw@ast.leeds.ac.uk), Tim Jeness (timj@jach.hawaii.edu),
 1638 and Doug Burke (burke@ifa.hawaii.edu).
 1639 
 1640 All rights reserved. There is no warranty. You are allowed
 1641 to redistribute this software / documentation under certain
 1642 conditions. For details, see the file COPYING in the PDL
 1643 distribution. If this file is separated from the PDL distribution,
 1644 the copyright notice should be included in the file.
 1645 
 1646 =cut
 1647 #line 1648 "Image2D.pm"
 1648 
 1649 
 1650 
 1651 
 1652 # Exit with OK status
 1653 
 1654 1;